Changeset 5246 for LMDZ6/trunk/libf/misc


Ignore:
Timestamp:
Oct 21, 2024, 2:58:45 PM (3 months ago)
Author:
abarral
Message:

Convert fixed-form to free-form sources .F -> .{f,F}90
(WIP: some .F remain, will be handled in subsequent commits)

Location:
LMDZ6/trunk/libf/misc
Files:
24 moved

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/misc/cbrt.f90

    r5245 r5246  
    1       FUNCTION cbrt(x)
    2       ! Return the cubic root for positive or negative x
    3       IMPLICIT NONE
     1FUNCTION cbrt(x)
     2  ! ! Return the cubic root for positive or negative x
     3  IMPLICIT NONE
    44
    5       REAL x,cbrt
     5  REAL :: x,cbrt
    66
    7       cbrt=sign(1.,x)*(abs(x)**(1./3.))
     7  cbrt=sign(1.,x)*(abs(x)**(1./3.))
    88
    9       RETURN
    10       END
     9  RETURN
     10END FUNCTION cbrt
    1111

  • LMDZ6/trunk/libf/misc/chfev.f90

    r5245 r5246  
    1 *DECK CHFEV
    2       SUBROUTINE CHFEV (X1, X2, F1, F2, D1, D2, NE, XE, FE, NEXT, IERR)
    3 C***BEGIN PROLOGUE  CHFEV
    4 C***PURPOSE  Evaluate a cubic polynomial given in Hermite form at an
    5 C            array of points.  While designed for use by PCHFE, it may
    6 C            be useful directly as an evaluator for a piecewise cubic
    7 C            Hermite function in applications, such as graphing, where
    8 C            the interval is known in advance.
    9 C***LIBRARY   SLATEC (PCHIP)
    10 C***CATEGORY  E3
    11 C***TYPE      SINGLE PRECISION (CHFEV-S, DCHFEV-D)
    12 C***KEYWORDS  CUBIC HERMITE EVALUATION, CUBIC POLYNOMIAL EVALUATION,
    13 C             PCHIP
    14 C***AUTHOR  Fritsch, F. N., (LLNL)
    15 C             Lawrence Livermore National Laboratory
    16 C             P.O. Box 808  (L-316)
    17 C             Livermore, CA  94550
    18 C             FTS 532-4275, (510) 422-4275
    19 C***DESCRIPTION
    20 C
    21 C          CHFEV:  Cubic Hermite Function EValuator
    22 C
    23 C    Evaluates the cubic polynomial determined by function values
    24 C    F1,F2 and derivatives D1,D2 on interval (X1,X2) at the points
    25 C    XE(J), J=1(1)NE.
    26 C
    27 C ----------------------------------------------------------------------
    28 C
    29 C  Calling sequence:
    30 C
    31 C        INTEGER  NE, NEXT(2), IERR
    32 C        REAL  X1, X2, F1, F2, D1, D2, XE(NE), FE(NE)
    33 C
    34 C        CALL  CHFEV (X1,X2, F1,F2, D1,D2, NE, XE, FE, NEXT, IERR)
    35 C
    36 C   Parameters:
    37 C
    38 C    X1,X2 -- (input) endpoints of interval of definition of cubic.
    39 C           (Error return if  X1.EQ.X2 .)
    40 C
    41 C    F1,F2 -- (input) values of function at X1 and X2, respectively.
    42 C
    43 C    D1,D2 -- (input) values of derivative at X1 and X2, respectively.
    44 C
    45 C    NE -- (input) number of evaluation points.  (Error return if
    46 C           NE.LT.1 .)
    47 C
    48 C    XE -- (input) real array of points at which the function is to be
    49 C           evaluated.  If any of the XE are outside the interval
    50 C           [X1,X2], a warning error is returned in NEXT.
    51 C
    52 C    FE -- (output) real array of values of the cubic function defined
    53 C           by  X1,X2, F1,F2, D1,D2  at the points  XE.
    54 C
    55 C    NEXT -- (output) integer array indicating number of extrapolation
    56 C           points:
    57 C            NEXT(1) = number of evaluation points to left of interval.
    58 C            NEXT(2) = number of evaluation points to right of interval.
    59 C
    60 C    IERR -- (output) error flag.
    61 C           Normal return:
    62 C              IERR = 0  (no errors).
    63 C           "Recoverable" errors:
    64 C              IERR = -1  if NE.LT.1 .
    65 C              IERR = -2  if X1.EQ.X2 .
    66 C                (The FE-array has not been changed in either case.)
    67 C
    68 C***REFERENCES  (NONE)
    69 C***ROUTINES CALLED  XERMSG
    70 C***REVISION HISTORY  (YYMMDD)
    71 C   811019  DATE WRITTEN
    72 C   820803  Minor cosmetic changes for release 1.
    73 C   890411  Added SAVE statements (Vers. 3.2).
    74 C   890531  Changed all specific intrinsics to generic.  (WRB)
    75 C   890703  Corrected category record.  (WRB)
    76 C   890703  REVISION DATE from Version 3.2
    77 C   891214  Prologue converted to Version 4.0 format.  (BAB)
    78 C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
    79 C***END PROLOGUE  CHFEV
    80 C  Programming notes:
    81 C
    82 C    To produce a double precision version, simply:
    83 C        a. Change CHFEV to DCHFEV wherever it occurs,
    84 C        b. Change the real declaration to double precision, and
    85 C        c. Change the constant ZERO to double precision.
    86 C
    87 C  DECLARE ARGUMENTS.
    88 C
    89       INTEGER NE, NEXT(2), IERR
    90       REAL X1, X2, F1, F2, D1, D2, XE(*), FE(*)
    91 C
    92 C  DECLARE LOCAL VARIABLES.
    93 C
    94       INTEGER I
    95       REAL C2, C3, DEL1, DEL2, DELTA, H, X, XMI, XMA, ZERO
    96       SAVE ZERO
    97       DATA  ZERO /0./
    98 C
    99 C  VALIDITY-CHECK ARGUMENTS.
    100 C
    101 C***FIRST EXECUTABLE STATEMENT  CHFEV
    102       IF (NE .LT. 1)  GO TO 5001
    103       H = X2 - X1
    104       IF (H .EQ. ZERO)  GO TO 5002
    105 C
    106 C  INITIALIZE.
    107 C
    108       IERR = 0
    109       NEXT(1) = 0
    110       NEXT(2) = 0
    111       XMI = MIN(ZERO, H)
    112       XMA = MAX(ZERO, H)
    113 C
    114 C  COMPUTE CUBIC COEFFICIENTS (EXPANDED ABOUT X1).
    115 C
    116       DELTA = (F2 - F1)/H
    117       DEL1 = (D1 - DELTA)/H
    118       DEL2 = (D2 - DELTA)/H
    119 C                                          (DELTA IS NO LONGER NEEDED.)
    120       C2 = -(DEL1+DEL1 + DEL2)
    121       C3 = (DEL1 + DEL2)/H
    122 C                              (H, DEL1 AND DEL2 ARE NO LONGER NEEDED.)
    123 C
    124 C  EVALUATION LOOP.
    125 C
    126       DO 500  I = 1, NE
    127          X = XE(I) - X1
    128          FE(I) = F1 + X*(D1 + X*(C2 + X*C3))
    129 C          COUNT EXTRAPOLATION POINTS.
    130          IF ( X.LT.XMI )  NEXT(1) = NEXT(1) + 1
    131          IF ( X.GT.XMA )  NEXT(2) = NEXT(2) + 1
    132 C        (NOTE REDUNDANCY--IF EITHER CONDITION IS TRUE, OTHER IS FALSE.)
    133   500 CONTINUE
    134 C
    135 C  NORMAL RETURN.
    136 C
    137       RETURN
    138 C
    139 C  ERROR RETURNS.
    140 C
    141  5001 CONTINUE
    142 C    NE.LT.1 RETURN.
    143       IERR = -1
    144       CALL XERMSG ('SLATEC', 'CHFEV',
    145      +   'NUMBER OF EVALUATION POINTS LESS THAN ONE', IERR, 1)
    146       RETURN
    147 C
    148  5002 CONTINUE
    149 C    X1.EQ.X2 RETURN.
    150       IERR = -2
    151       CALL XERMSG ('SLATEC', 'CHFEV', 'INTERVAL ENDPOINTS EQUAL', IERR,
    152      +   1)
    153       RETURN
    154 C------------- LAST LINE OF CHFEV FOLLOWS ------------------------------
    155       END
     1!DECK CHFEV
     2SUBROUTINE CHFEV (X1, X2, F1, F2, D1, D2, NE, XE, FE, NEXT, IERR)
     3  !***BEGIN PROLOGUE  CHFEV
     4  !***PURPOSE  Evaluate a cubic polynomial given in Hermite form at an
     5         ! array of points.  While designed for use by PCHFE, it may
     6         ! be useful directly as an evaluator for a piecewise cubic
     7         ! Hermite function in applications, such as graphing, where
     8         ! the interval is known in advance.
     9  !***LIBRARY   SLATEC (PCHIP)
     10  !***CATEGORY  E3
     11  !***TYPE      SINGLE PRECISION (CHFEV-S, DCHFEV-D)
     12  !***KEYWORDS  CUBIC HERMITE EVALUATION, CUBIC POLYNOMIAL EVALUATION,
     13         !  PCHIP
     14  !***AUTHOR  Fritsch, F. N., (LLNL)
     15         !  Lawrence Livermore National Laboratory
     16         !  P.O. Box 808  (L-316)
     17         !  Livermore, CA  94550
     18         !  FTS 532-4275, (510) 422-4275
     19  !***DESCRIPTION
     20  !
     21  !      CHFEV:  Cubic Hermite Function EValuator
     22  !
     23  ! Evaluates the cubic polynomial determined by function values
     24  ! F1,F2 and derivatives D1,D2 on interval (X1,X2) at the points
     25  ! XE(J), J=1(1)NE.
     26  !
     27  ! ----------------------------------------------------------------------
     28  !
     29  !  Calling sequence:
     30  !
     31  !    INTEGER  NE, NEXT(2), IERR
     32  !    REAL  X1, X2, F1, F2, D1, D2, XE(NE), FE(NE)
     33  !
     34  !    CALL  CHFEV (X1,X2, F1,F2, D1,D2, NE, XE, FE, NEXT, IERR)
     35  !
     36  !   Parameters:
     37  !
     38  ! X1,X2 -- (input) endpoints of interval of definition of cubic.
     39  !       (Error return if  X1.EQ.X2 .)
     40  !
     41  ! F1,F2 -- (input) values of function at X1 and X2, respectively.
     42  !
     43  ! D1,D2 -- (input) values of derivative at X1 and X2, respectively.
     44  !
     45  ! NE -- (input) number of evaluation points.  (Error return if
     46  !       NE.LT.1 .)
     47  !
     48  ! XE -- (input) real array of points at which the function is to be
     49  !       evaluated.  If any of the XE are outside the interval
     50  !       [X1,X2], a warning error is returned in NEXT.
     51  !
     52  ! FE -- (output) real array of values of the cubic function defined
     53  !       by  X1,X2, F1,F2, D1,D2  at the points  XE.
     54  !
     55  ! NEXT -- (output) integer array indicating number of extrapolation
     56  !       points:
     57  !        NEXT(1) = number of evaluation points to left of interval.
     58  !        NEXT(2) = number of evaluation points to right of interval.
     59  !
     60  ! IERR -- (output) error flag.
     61  !       Normal return:
     62  !          IERR = 0  (no errors).
     63  !       "Recoverable" errors:
     64  !          IERR = -1  if NE.LT.1 .
     65  !          IERR = -2  if X1.EQ.X2 .
     66  !            (The FE-array has not been changed in either case.)
     67  !
     68  !***REFERENCES  (NONE)
     69  !***ROUTINES CALLED  XERMSG
     70  !***REVISION HISTORY  (YYMMDD)
     71  !   811019  DATE WRITTEN
     72  !   820803  Minor cosmetic changes for release 1.
     73  !   890411  Added SAVE statements (Vers. 3.2).
     74  !   890531  Changed all specific intrinsics to generic.  (WRB)
     75  !   890703  Corrected category record.  (WRB)
     76  !   890703  REVISION DATE from Version 3.2
     77  !   891214  Prologue converted to Version 4.0 format.  (BAB)
     78  !   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
     79  !***END PROLOGUE  CHFEV
     80  !  Programming notes:
     81  !
     82  ! To produce a double precision version, simply:
     83  !    a. Change CHFEV to DCHFEV wherever it occurs,
     84  !    b. Change the real declaration to double precision, and
     85  !    c. Change the constant ZERO to double precision.
     86  !
     87  !  DECLARE ARGUMENTS.
     88  !
     89  INTEGER :: NE, NEXT(2), IERR
     90  REAL :: X1, X2, F1, F2, D1, D2, XE(*), FE(*)
     91  !
     92  !  DECLARE LOCAL VARIABLES.
     93  !
     94  INTEGER :: I
     95  REAL :: C2, C3, DEL1, DEL2, DELTA, H, X, XMI, XMA, ZERO
     96  SAVE ZERO
     97  DATA  ZERO /0./
     98  !
     99  !  VALIDITY-CHECK ARGUMENTS.
     100  !
     101  !***FIRST EXECUTABLE STATEMENT  CHFEV
     102  IF (NE .LT. 1)  GO TO 5001
     103  H = X2 - X1
     104  IF (H .EQ. ZERO)  GO TO 5002
     105  !
     106  !  INITIALIZE.
     107  !
     108  IERR = 0
     109  NEXT(1) = 0
     110  NEXT(2) = 0
     111  XMI = MIN(ZERO, H)
     112  XMA = MAX(ZERO, H)
     113  !
     114  !  COMPUTE CUBIC COEFFICIENTS (EXPANDED ABOUT X1).
     115  !
     116  DELTA = (F2 - F1)/H
     117  DEL1 = (D1 - DELTA)/H
     118  DEL2 = (D2 - DELTA)/H
     119                                        ! (DELTA IS NO LONGER NEEDED.)
     120  C2 = -(DEL1+DEL1 + DEL2)
     121  C3 = (DEL1 + DEL2)/H
     122                            ! (H, DEL1 AND DEL2 ARE NO LONGER NEEDED.)
     123  !
     124  !  EVALUATION LOOP.
     125  !
     126  DO  I = 1, NE
     127     X = XE(I) - X1
     128     FE(I) = F1 + X*(D1 + X*(C2 + X*C3))
     129       ! COUNT EXTRAPOLATION POINTS.
     130     IF ( X.LT.XMI )  NEXT(1) = NEXT(1) + 1
     131     IF ( X.GT.XMA )  NEXT(2) = NEXT(2) + 1
     132     ! (NOTE REDUNDANCY--IF EITHER CONDITION IS TRUE, OTHER IS FALSE.)
     133  END DO
     134  !
     135  !  NORMAL RETURN.
     136  !
     137  RETURN
     138  !
     139  !  ERROR RETURNS.
     140  !
     141 5001   CONTINUE
     142  ! NE.LT.1 RETURN.
     143  IERR = -1
     144  CALL XERMSG ('SLATEC', 'CHFEV', &
     145        'NUMBER OF EVALUATION POINTS LESS THAN ONE', IERR, 1)
     146  RETURN
     147  !
     148 5002   CONTINUE
     149  ! X1.EQ.X2 RETURN.
     150  IERR = -2
     151  CALL XERMSG ('SLATEC', 'CHFEV', 'INTERVAL ENDPOINTS EQUAL', IERR, &
     152        1)
     153  RETURN
     154  !------------- LAST LINE OF CHFEV FOLLOWS ------------------------------
     155END SUBROUTINE CHFEV

  • LMDZ6/trunk/libf/misc/cray.F90

    r5245 r5246  
    33!
    44#ifdef CRAY
    5       SUBROUTINE riencray
    6       END
     5SUBROUTINE riencray
     6END SUBROUTINE riencray
    77#else
    8       subroutine scopy(n,sx,incx,sy,incy)
    9 c
    10       IMPLICIT NONE
    11 c
    12       integer n,incx,incy,ix,iy,i
    13       real sx((n-1)*incx+1),sy((n-1)*incy+1)
    14 c
    15       iy=1
    16       ix=1
    17       do 10 i=1,n
    18          sy(iy)=sx(ix)
    19          ix=ix+incx
    20          iy=iy+incy
    21 10    continue
    22 c
    23       return
    24       end
     8subroutine scopy(n,sx,incx,sy,incy)
     9  !
     10  IMPLICIT NONE
     11  !
     12  integer :: n,incx,incy,ix,iy,i
     13  real :: sx((n-1)*incx+1),sy((n-1)*incy+1)
     14  !
     15  iy=1
     16  ix=1
     17  do i=1,n
     18     sy(iy)=sx(ix)
     19     ix=ix+incx
     20     iy=iy+incy
     21  end do
     22  !
     23  return
     24end subroutine scopy
    2525
    26       function ssum(n,sx,incx)
    27 c
    28       IMPLICIT NONE
    29 c
    30       integer n,incx,i,ix
    31       real ssum,sx((n-1)*incx+1)
    32 c
    33       ssum=0.
    34       ix=1
    35       do 10 i=1,n
    36          ssum=ssum+sx(ix)
    37          ix=ix+incx
    38 10    continue
    39 c
    40       return
    41       end
     26function ssum(n,sx,incx)
     27  !
     28  IMPLICIT NONE
     29  !
     30  integer :: n,incx,i,ix
     31  real :: ssum,sx((n-1)*incx+1)
     32  !
     33  ssum=0.
     34  ix=1
     35  do i=1,n
     36     ssum=ssum+sx(ix)
     37     ix=ix+incx
     38  end do
     39  !
     40  return
     41end function ssum
    4242#endif

  • LMDZ6/trunk/libf/misc/fdump.f90

    r5245 r5246  
    1 *DECK FDUMP
    2       SUBROUTINE FDUMP
    3 C***BEGIN PROLOGUE  FDUMP
    4 C***PURPOSE  Symbolic dump (should be locally written).
    5 C***LIBRARY   SLATEC (XERROR)
    6 C***CATEGORY  R3
    7 C***TYPE      ALL (FDUMP-A)
    8 C***KEYWORDS  ERROR, XERMSG
    9 C***AUTHOR  Jones, R. E., (SNLA)
    10 C***DESCRIPTION
    11 C
    12 C        ***Note*** Machine Dependent Routine
    13 C        FDUMP is intended to be replaced by a locally written
    14 C        version which produces a symbolic dump.  Failing this,
    15 C        it should be replaced by a version which prints the
    16 C        subprogram nesting list.  Note that this dump must be
    17 C        printed on each of up to five files, as indicated by the
    18 C        XGETUA routine.  See XSETUA and XGETUA for details.
    19 C
    20 C    Written by Ron Jones, with SLATEC Common Math Library Subcommittee
    21 C
    22 C***REFERENCES  (NONE)
    23 C***ROUTINES CALLED  (NONE)
    24 C***REVISION HISTORY  (YYMMDD)
    25 C   790801  DATE WRITTEN
    26 C   861211  REVISION DATE from Version 3.2
    27 C   891214  Prologue converted to Version 4.0 format.  (BAB)
    28 C***END PROLOGUE  FDUMP
    29 C***FIRST EXECUTABLE STATEMENT  FDUMP
    30       RETURN
    31       END
     1!DECK FDUMP
     2SUBROUTINE FDUMP
     3  !***BEGIN PROLOGUE  FDUMP
     4  !***PURPOSE  Symbolic dump (should be locally written).
     5  !***LIBRARY   SLATEC (XERROR)
     6  !***CATEGORY  R3
     7  !***TYPE      ALL (FDUMP-A)
     8  !***KEYWORDS  ERROR, XERMSG
     9  !***AUTHOR  Jones, R. E., (SNLA)
     10  !***DESCRIPTION
     11  !
     12  !    ***Note*** Machine Dependent Routine
     13  !    FDUMP is intended to be replaced by a locally written
     14  !    version which produces a symbolic dump.  Failing this,
     15  !    it should be replaced by a version which prints the
     16  !    subprogram nesting list.  Note that this dump must be
     17  !    printed on each of up to five files, as indicated by the
     18  !    XGETUA routine.  See XSETUA and XGETUA for details.
     19  !
     20  ! Written by Ron Jones, with SLATEC Common Math Library Subcommittee
     21  !
     22  !***REFERENCES  (NONE)
     23  !***ROUTINES CALLED  (NONE)
     24  !***REVISION HISTORY  (YYMMDD)
     25  !   790801  DATE WRITTEN
     26  !   861211  REVISION DATE from Version 3.2
     27  !   891214  Prologue converted to Version 4.0 format.  (BAB)
     28  !***END PROLOGUE  FDUMP
     29  !***FIRST EXECUTABLE STATEMENT  FDUMP
     30  RETURN
     31END SUBROUTINE FDUMP

  • LMDZ6/trunk/libf/misc/formcoord.f90

    r5245 r5246  
    22! $Header$
    33!
    4       subroutine formcoord(unit,n,x,a,rev,text)
    5       implicit none
    6       integer n,unit,ndec
    7       logical rev
    8       real x(n),a
    9       character*4 text
     4subroutine formcoord(unit,n,x,a,rev,text)
     5  implicit none
     6  integer :: n,unit,ndec
     7  logical :: rev
     8  real :: x(n),a
     9  character(len=4) :: text
    1010
    11       integer i,id,i1,i2,in
    12       real dx,dxmin
     11  integer :: i,id,i1,i2,in
     12  real :: dx,dxmin
    1313
    14       if(rev) then
    15          id=-1
    16          i1=n
    17          i2=n-1
    18          in=1
    19          write(unit,3000) text(1:1)
    20       else
    21          id=1
    22          i1=1
    23          i2=2
    24          in=n
    25       endif
     14  if(rev) then
     15     id=-1
     16     i1=n
     17     i2=n-1
     18     in=1
     19     write(unit,3000) text(1:1)
     20  else
     21     id=1
     22     i1=1
     23     i2=2
     24     in=n
     25  endif
    2626
    27       if (n.lt.2) then
    28          ndec=1
    29          write(unit,1000) text,n,x(1)*a
    30       else
    31          dxmin=abs(x(2)-x(1))
    32          do i=2,n-1
    33             dx=abs(x(i+1)-x(i))
    34             if (dx.lt.dxmin) dxmin=dx
    35          enddo
     27  if (n.lt.2) then
     28     ndec=1
     29     write(unit,1000) text,n,x(1)*a
     30  else
     31     dxmin=abs(x(2)-x(1))
     32     do i=2,n-1
     33        dx=abs(x(i+1)-x(i))
     34        if (dx.lt.dxmin) dxmin=dx
     35     enddo
    3636
    37          ndec=-log10(dxmin)+2
    38          if(mod(n,6).eq.1) then
    39             write(unit,1000) text,n,x(i1)*a
    40             write(unit,2000) (x(i)*a,i=i2,in,id)
    41          else
    42             write(unit,1000) text,n
    43             write(unit,2000) (x(i)*a,i=i1,in,id)
    44          endif
    45       endif
     37     ndec=-log10(dxmin)+2
     38     if(mod(n,6).eq.1) then
     39        write(unit,1000) text,n,x(i1)*a
     40        write(unit,2000) (x(i)*a,i=i2,in,id)
     41     else
     42        write(unit,1000) text,n
     43        write(unit,2000) (x(i)*a,i=i1,in,id)
     44     endif
     45  endif
    4646
    47 1000  format(a4,2x,i4,' LEVELS',43x,f12.2)
    48 2000  format(6f12.2)
    49 c1000  format(a4,2x,i4,' LEVELS',43x,f12.<ndec>)
    50 c2000  format(6f12.<ndec>)
    51 3000  format('FORMAT ',a1,'REV')
    52       return
     471000   format(a4,2x,i4,' LEVELS',43x,f12.2)
     482000   format(6f12.2)
     49  !1000  format(a4,2x,i4,' LEVELS',43x,f12.<ndec>)
     50  !2000  format(6f12.<ndec>)
     513000   format('FORMAT ',a1,'REV')
     52  return
    5353
    54       end
     54end subroutine formcoord

  • LMDZ6/trunk/libf/misc/i1mach.f90

    r5245 r5246  
    1 *DECK I1MACH
    2       INTEGER FUNCTION I1MACH (I)
    3       IMPLICIT NONE
    4 C***BEGIN PROLOGUE  I1MACH
    5 C***PURPOSE  Return integer machine dependent constants.
    6 C***LIBRARY   SLATEC
    7 C***CATEGORY  R1
    8 C***TYPE      INTEGER (I1MACH-I)
    9 C***KEYWORDS  MACHINE CONSTANTS
    10 C***AUTHOR  Fox, P. A., (Bell Labs)
    11 C          Hall, A. D., (Bell Labs)
    12 C          Schryer, N. L., (Bell Labs)
    13 C***DESCRIPTION
    14 C
    15 C   I1MACH can be used to obtain machine-dependent parameters for the
    16 C   local machine environment.  It is a function subprogram with one
    17 C   (input) argument and can be referenced as follows:
    18 C
    19 C        K = I1MACH(I)
    20 C
    21 C   where I=1,...,16.  The (output) value of K above is determined by
    22 C   the (input) value of I.  The results for various values of I are
    23 C   discussed below.
    24 C
    25 C   I/O unit numbers:
    26 C    I1MACH( 1) = the standard input unit.
    27 C    I1MACH( 2) = the standard output unit.
    28 C    I1MACH( 3) = the standard punch unit.
    29 C    I1MACH( 4) = the standard error message unit.
    30 C
    31 C   Words:
    32 C    I1MACH( 5) = the number of bits per integer storage unit.
    33 C    I1MACH( 6) = the number of characters per integer storage unit.
    34 C
    35 C   Integers:
    36 C    assume integers are represented in the S-digit, base-A form
    37 C
    38 C                sign ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) )
    39 C
    40 C                where 0 .LE. X(I) .LT. A for I=0,...,S-1.
    41 C    I1MACH( 7) = A, the base.
    42 C    I1MACH( 8) = S, the number of base-A digits.
    43 C    I1MACH( 9) = A**S - 1, the largest magnitude.
    44 C
    45 C   Floating-Point Numbers:
    46 C    Assume floating-point numbers are represented in the T-digit,
    47 C    base-B form
    48 C                sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
    49 C
    50 C                where 0 .LE. X(I) .LT. B for I=1,...,T,
    51 C                0 .LT. X(1), and EMIN .LE. E .LE. EMAX.
    52 C    I1MACH(10) = B, the base.
    53 C
    54 C   Single-Precision:
    55 C    I1MACH(11) = T, the number of base-B digits.
    56 C    I1MACH(12) = EMIN, the smallest exponent E.
    57 C    I1MACH(13) = EMAX, the largest exponent E.
    58 C
    59 C   Double-Precision:
    60 C    I1MACH(14) = T, the number of base-B digits.
    61 C    I1MACH(15) = EMIN, the smallest exponent E.
    62 C    I1MACH(16) = EMAX, the largest exponent E.
    63 C
    64 C   To alter this function for a particular environment, the desired
    65 C   set of DATA statements should be activated by removing the C from
    66 C   column 1.  Also, the values of I1MACH(1) - I1MACH(4) should be
    67 C   checked for consistency with the local operating system.
    68 C
    69 C***REFERENCES  P. A. Fox, A. D. Hall and N. L. Schryer, Framework for
    70 C                 a portable library, ACM Transactions on Mathematical
    71 C                 Software 4, 2 (June 1978), pp. 177-188.
    72 C***ROUTINES CALLED  (NONE)
    73 C***REVISION HISTORY  (YYMMDD)
    74 C   750101  DATE WRITTEN
    75 C   891012  Added VAX G-floating constants.  (WRB)
    76 C   891012  REVISION DATE from Version 3.2
    77 C   891214  Prologue converted to Version 4.0 format.  (BAB)
    78 C   900618  Added DEC RISC constants.  (WRB)
    79 C   900723  Added IBM RS 6000 constants.  (WRB)
    80 C   901009  Correct I1MACH(7) for IBM Mainframes. Should be 2 not 16.
    81 C           (RWC)
    82 C   910710  Added HP 730 constants.  (SMR)
    83 C   911114  Added Convex IEEE constants.  (WRB)
    84 C   920121  Added SUN -r8 compiler option constants.  (WRB)
    85 C   920229  Added Touchstone Delta i860 constants.  (WRB)
    86 C   920501  Reformatted the REFERENCES section.  (WRB)
    87 C   920625  Added Convex -p8 and -pd8 compiler option constants.
    88 C           (BKS, WRB)
    89 C   930201  Added DEC Alpha and SGI constants.  (RWC and WRB)
    90 C   930618  Corrected I1MACH(5) for Convex -p8 and -pd8 compiler
    91 C           options.  (DWL, RWC and WRB).
    92 C   100623  Use Fortran 95 intrinsic functions (Lionel GUEZ)
    93 C***END PROLOGUE  I1MACH
    94 C
    95       INTEGER IMACH(16),OUTPUT
    96       SAVE IMACH
    97       EQUIVALENCE (IMACH(4),OUTPUT)
    98       INTEGER I
    99 C***FIRST EXECUTABLE STATEMENT  I1MACH
    100       IMACH( 1) =         5
    101       IMACH( 2) =         6
    102       IMACH( 3) =         6
    103       IMACH( 4) =         6
    104       IMACH( 5) =        bit_size(0)
    105       IMACH( 6) =         IMACH( 5) / 8
    106       IMACH( 7) =         radix(0)
    107       IMACH( 8) =        digits(0)
    108       IMACH( 9) =     huge(0)
    109       IMACH(10) =         radix(0.)
    110       IMACH(11) =        digits(0.)
    111       IMACH(12) =      minexponent(0.)
    112       IMACH(13) =       maxexponent(0.)
    113       IMACH(14) =        digits(0d0)
    114       IMACH(15) =      minexponent(0d0)
    115       IMACH(16) =       maxexponent(0d0)
    116       IF (I .LT. 1  .OR.  I .GT. 16) GO TO 10
    117 C
    118       I1MACH = IMACH(I)
    119       RETURN
    120 C
    121    10 CONTINUE
    122       WRITE (UNIT = OUTPUT, FMT = 9000)
    123  9000 FORMAT ('1ERROR    1 IN I1MACH - I OUT OF BOUNDS')
    124 C
    125 C    CALL FDUMP
    126 C
    127       STOP
    128       END
     1!DECK I1MACH
     2INTEGER FUNCTION I1MACH (I)
     3  IMPLICIT NONE
     4  !***BEGIN PROLOGUE  I1MACH
     5  !***PURPOSE  Return integer machine dependent constants.
     6  !***LIBRARY   SLATEC
     7  !***CATEGORY  R1
     8  !***TYPE      INTEGER (I1MACH-I)
     9  !***KEYWORDS  MACHINE CONSTANTS
     10  !***AUTHOR  Fox, P. A., (Bell Labs)
     11        ! Hall, A. D., (Bell Labs)
     12        ! Schryer, N. L., (Bell Labs)
     13  !***DESCRIPTION
     14  !
     15  !   I1MACH can be used to obtain machine-dependent parameters for the
     16  !   local machine environment.  It is a function subprogram with one
     17  !   (input) argument and can be referenced as follows:
     18  !
     19  !    K = I1MACH(I)
     20  !
     21  !   where I=1,...,16.  The (output) value of K above is determined by
     22  !   the (input) value of I.  The results for various values of I are
     23  !   discussed below.
     24  !
     25  !   I/O unit numbers:
     26  ! I1MACH( 1) = the standard input unit.
     27  ! I1MACH( 2) = the standard output unit.
     28  ! I1MACH( 3) = the standard punch unit.
     29  ! I1MACH( 4) = the standard error message unit.
     30  !
     31  !   Words:
     32  ! I1MACH( 5) = the number of bits per integer storage unit.
     33  ! I1MACH( 6) = the number of characters per integer storage unit.
     34  !
     35  !   Integers:
     36  ! assume integers are represented in the S-digit, base-A form
     37  !
     38  !            sign ( X(S-1)*A**(S-1) + ... + X(1)*A + X(0) )
     39  !
     40  !            where 0 .LE. X(I) .LT. A for I=0,...,S-1.
     41  ! I1MACH( 7) = A, the base.
     42  ! I1MACH( 8) = S, the number of base-A digits.
     43  ! I1MACH( 9) = A**S - 1, the largest magnitude.
     44  !
     45  !   Floating-Point Numbers:
     46  ! Assume floating-point numbers are represented in the T-digit,
     47  ! base-B form
     48  !            sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
     49  !
     50  !            where 0 .LE. X(I) .LT. B for I=1,...,T,
     51  !            0 .LT. X(1), and EMIN .LE. E .LE. EMAX.
     52  ! I1MACH(10) = B, the base.
     53  !
     54  !   Single-Precision:
     55  ! I1MACH(11) = T, the number of base-B digits.
     56  ! I1MACH(12) = EMIN, the smallest exponent E.
     57  ! I1MACH(13) = EMAX, the largest exponent E.
     58  !
     59  !   Double-Precision:
     60  ! I1MACH(14) = T, the number of base-B digits.
     61  ! I1MACH(15) = EMIN, the smallest exponent E.
     62  ! I1MACH(16) = EMAX, the largest exponent E.
     63  !
     64  !   To alter this function for a particular environment, the desired
     65  !   set of DATA statements should be activated by removing the C from
     66  !   column 1.  Also, the values of I1MACH(1) - I1MACH(4) should be
     67  !   checked for consistency with the local operating system.
     68  !
     69  !***REFERENCES  P. A. Fox, A. D. Hall and N. L. Schryer, Framework for
     70  !             a portable library, ACM Transactions on Mathematical
     71  !             Software 4, 2 (June 1978), pp. 177-188.
     72  !***ROUTINES CALLED  (NONE)
     73  !***REVISION HISTORY  (YYMMDD)
     74  !   750101  DATE WRITTEN
     75  !   891012  Added VAX G-floating constants.  (WRB)
     76  !   891012  REVISION DATE from Version 3.2
     77  !   891214  Prologue converted to Version 4.0 format.  (BAB)
     78  !   900618  Added DEC RISC constants.  (WRB)
     79  !   900723  Added IBM RS 6000 constants.  (WRB)
     80  !   901009  Correct I1MACH(7) for IBM Mainframes. Should be 2 not 16.
     81  !       (RWC)
     82  !   910710  Added HP 730 constants.  (SMR)
     83  !   911114  Added Convex IEEE constants.  (WRB)
     84  !   920121  Added SUN -r8 compiler option constants.  (WRB)
     85  !   920229  Added Touchstone Delta i860 constants.  (WRB)
     86  !   920501  Reformatted the REFERENCES section.  (WRB)
     87  !   920625  Added Convex -p8 and -pd8 compiler option constants.
     88  !       (BKS, WRB)
     89  !   930201  Added DEC Alpha and SGI constants.  (RWC and WRB)
     90  !   930618  Corrected I1MACH(5) for Convex -p8 and -pd8 compiler
     91  !       options.  (DWL, RWC and WRB).
     92  !   100623  Use Fortran 95 intrinsic functions (Lionel GUEZ)
     93  !***END PROLOGUE  I1MACH
     94  !
     95  INTEGER :: IMACH(16),OUTPUT
     96  SAVE IMACH
     97  EQUIVALENCE (IMACH(4),OUTPUT)
     98  INTEGER :: I
     99  !***FIRST EXECUTABLE STATEMENT  I1MACH
     100  IMACH( 1) =         5
     101  IMACH( 2) =         6
     102  IMACH( 3) =         6
     103  IMACH( 4) =         6
     104  IMACH( 5) =        bit_size(0)
     105  IMACH( 6) =         IMACH( 5) / 8
     106  IMACH( 7) =         radix(0)
     107  IMACH( 8) =        digits(0)
     108  IMACH( 9) =     huge(0)
     109  IMACH(10) =         radix(0.)
     110  IMACH(11) =        digits(0.)
     111  IMACH(12) =      minexponent(0.)
     112  IMACH(13) =       maxexponent(0.)
     113  IMACH(14) =        digits(0d0)
     114  IMACH(15) =      minexponent(0d0)
     115  IMACH(16) =       maxexponent(0d0)
     116  IF (I .LT. 1  .OR.  I .GT. 16) GO TO 10
     117  !
     118  I1MACH = IMACH(I)
     119  RETURN
     120  !
     121   10   CONTINUE
     122  WRITE (UNIT = OUTPUT, FMT = 9000)
     123 9000   FORMAT ('1ERROR    1 IN I1MACH - I OUT OF BOUNDS')
     124  !
     125  ! CALL FDUMP
     126  !
     127  STOP
     128END FUNCTION I1MACH

  • LMDZ6/trunk/libf/misc/ismax.f90

    r5245 r5246  
    22! $Header$
    33!
    4       function ismax(n,sx,incx)
    5 c
    6       IMPLICIT NONE
    7 c
    8       INTEGER n,i,incx,ismax,ix
    9       real sx((n-1)*incx+1),sxmax
    10 c
    11       ix=1
    12       ismax=1
    13       sxmax=sx(1)
    14       do 10 i=1,n-1
    15        ix=ix+incx
    16        if(sx(ix).gt.sxmax) then
    17          sxmax=sx(ix)
    18          ismax=i+1
    19        endif
    20 10    continue
    21 c
    22       return
    23       end
     4function ismax(n,sx,incx)
     5  !
     6  IMPLICIT NONE
     7  !
     8  INTEGER :: n,i,incx,ismax,ix
     9  real :: sx((n-1)*incx+1),sxmax
     10  !
     11  ix=1
     12  ismax=1
     13  sxmax=sx(1)
     14  do i=1,n-1
     15   ix=ix+incx
     16   if(sx(ix).gt.sxmax) then
     17     sxmax=sx(ix)
     18     ismax=i+1
     19   endif
     20  end do
     21  !
     22  return
     23end function ismax
    2424

  • LMDZ6/trunk/libf/misc/ismin.f90

    r5245 r5246  
    22! $Header$
    33!
    4       FUNCTION ismin(n,sx,incx)
    5 c
    6       IMPLICIT NONE
    7 c
    8       integer n,i,incx,ismin,ix
    9       real sx((n-1)*incx+1),sxmin
    10 c
    11       ix=1
    12       ismin=1
    13       sxmin=sx(1)
    14       DO i=1,n-1
    15          ix=ix+incx
    16          if(sx(ix).lt.sxmin) then
    17              sxmin=sx(ix)
    18              ismin=i+1
    19          endif
    20       ENDDO
    21 c
    22       return
    23       end
    24 C
     4FUNCTION ismin(n,sx,incx)
     5  !
     6  IMPLICIT NONE
     7  !
     8  integer :: n,i,incx,ismin,ix
     9  real :: sx((n-1)*incx+1),sxmin
     10  !
     11  ix=1
     12  ismin=1
     13  sxmin=sx(1)
     14  DO i=1,n-1
     15     ix=ix+incx
     16     if(sx(ix).lt.sxmin) then
     17         sxmin=sx(ix)
     18         ismin=i+1
     19     endif
     20  ENDDO
     21  !
     22  return
     23end function ismin
     24!

  • LMDZ6/trunk/libf/misc/j4save.f90

    r5245 r5246  
    1 *DECK J4SAVE
    2       FUNCTION J4SAVE (IWHICH, IVALUE, ISET)
    3       IMPLICIT NONE
    4 C***BEGIN PROLOGUE  J4SAVE
    5 C***SUBSIDIARY
    6 C***PURPOSE  Save or recall global variables needed by error
    7 C            handling routines.
    8 C***LIBRARY   SLATEC (XERROR)
    9 C***TYPE      INTEGER (J4SAVE-I)
    10 C***KEYWORDS  ERROR MESSAGES, ERROR NUMBER, RECALL, SAVE, XERROR
    11 C***AUTHOR  Jones, R. E., (SNLA)
    12 C***DESCRIPTION
    13 C
    14 C    Abstract
    15 C        J4SAVE saves and recalls several global variables needed
    16 C        by the library error handling routines.
    17 C
    18 C    Description of Parameters
    19 C      --Input--
    20 C        IWHICH - Index of item desired.
    21 C                = 1 Refers to current error number.
    22 C                = 2 Refers to current error control flag.
    23 C                = 3 Refers to current unit number to which error
    24 C                    messages are to be sent.  (0 means use standard.)
    25 C                = 4 Refers to the maximum number of times any
    26 C                     message is to be printed (as set by XERMAX).
    27 C                = 5 Refers to the total number of units to which
    28 C                     each error message is to be written.
    29 C                = 6 Refers to the 2nd unit for error messages
    30 C                = 7 Refers to the 3rd unit for error messages
    31 C                = 8 Refers to the 4th unit for error messages
    32 C                = 9 Refers to the 5th unit for error messages
    33 C        IVALUE - The value to be set for the IWHICH-th parameter,
    34 C                 if ISET is .TRUE. .
    35 C        ISET   - If ISET=.TRUE., the IWHICH-th parameter will BE
    36 C                 given the value, IVALUE.  If ISET=.FALSE., the
    37 C                 IWHICH-th parameter will be unchanged, and IVALUE
    38 C                 is a dummy parameter.
    39 C      --Output--
    40 C        The (old) value of the IWHICH-th parameter will be returned
    41 C        in the function value, J4SAVE.
    42 C
    43 C***SEE ALSO  XERMSG
    44 C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
    45 C                 Error-handling Package, SAND82-0800, Sandia
    46 C                 Laboratories, 1982.
    47 C***ROUTINES CALLED  (NONE)
    48 C***REVISION HISTORY  (YYMMDD)
    49 C   790801  DATE WRITTEN
    50 C   891214  Prologue converted to Version 4.0 format.  (BAB)
    51 C   900205  Minor modifications to prologue.  (WRB)
    52 C   900402  Added TYPE section.  (WRB)
    53 C   910411  Added KEYWORDS section.  (WRB)
    54 C   920501  Reformatted the REFERENCES section.  (WRB)
    55 C***END PROLOGUE  J4SAVE
    56       LOGICAL ISET
    57       INTEGER IPARAM(9)
    58       SAVE IPARAM
    59       DATA IPARAM(1),IPARAM(2),IPARAM(3),IPARAM(4)/0,2,0,10/
    60       DATA IPARAM(5)/1/
    61       DATA IPARAM(6),IPARAM(7),IPARAM(8),IPARAM(9)/0,0,0,0/
    62       INTEGER J4SAVE,IWHICH,IVALUE
    63 C***FIRST EXECUTABLE STATEMENT  J4SAVE
    64       J4SAVE = IPARAM(IWHICH)
    65       IF (ISET) IPARAM(IWHICH) = IVALUE
    66       RETURN
    67       END
     1!DECK J4SAVE
     2FUNCTION J4SAVE (IWHICH, IVALUE, ISET)
     3  IMPLICIT NONE
     4  !***BEGIN PROLOGUE  J4SAVE
     5  !***SUBSIDIARY
     6  !***PURPOSE  Save or recall global variables needed by error
     7         ! handling routines.
     8  !***LIBRARY   SLATEC (XERROR)
     9  !***TYPE      INTEGER (J4SAVE-I)
     10  !***KEYWORDS  ERROR MESSAGES, ERROR NUMBER, RECALL, SAVE, XERROR
     11  !***AUTHOR  Jones, R. E., (SNLA)
     12  !***DESCRIPTION
     13  !
     14  ! Abstract
     15  !    J4SAVE saves and recalls several global variables needed
     16  !    by the library error handling routines.
     17  !
     18  ! Description of Parameters
     19  !  --Input--
     20  !    IWHICH - Index of item desired.
     21  !            = 1 Refers to current error number.
     22  !            = 2 Refers to current error control flag.
     23  !            = 3 Refers to current unit number to which error
     24  !                messages are to be sent.  (0 means use standard.)
     25  !            = 4 Refers to the maximum number of times any
     26  !                 message is to be printed (as set by XERMAX).
     27  !            = 5 Refers to the total number of units to which
     28  !                 each error message is to be written.
     29  !            = 6 Refers to the 2nd unit for error messages
     30  !            = 7 Refers to the 3rd unit for error messages
     31  !            = 8 Refers to the 4th unit for error messages
     32  !            = 9 Refers to the 5th unit for error messages
     33  !    IVALUE - The value to be set for the IWHICH-th parameter,
     34  !             if ISET is .TRUE. .
     35  !    ISET   - If ISET=.TRUE., the IWHICH-th parameter will BE
     36  !             given the value, IVALUE.  If ISET=.FALSE., the
     37  !             IWHICH-th parameter will be unchanged, and IVALUE
     38  !             is a dummy parameter.
     39  !  --Output--
     40  !    The (old) value of the IWHICH-th parameter will be returned
     41  !    in the function value, J4SAVE.
     42  !
     43  !***SEE ALSO  XERMSG
     44  !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
     45  !             Error-handling Package, SAND82-0800, Sandia
     46  !             Laboratories, 1982.
     47  !***ROUTINES CALLED  (NONE)
     48  !***REVISION HISTORY  (YYMMDD)
     49  !   790801  DATE WRITTEN
     50  !   891214  Prologue converted to Version 4.0 format.  (BAB)
     51  !   900205  Minor modifications to prologue.  (WRB)
     52  !   900402  Added TYPE section.  (WRB)
     53  !   910411  Added KEYWORDS section.  (WRB)
     54  !   920501  Reformatted the REFERENCES section.  (WRB)
     55  !***END PROLOGUE  J4SAVE
     56  LOGICAL :: ISET
     57  INTEGER :: IPARAM(9)
     58  SAVE IPARAM
     59  DATA IPARAM(1),IPARAM(2),IPARAM(3),IPARAM(4)/0,2,0,10/
     60  DATA IPARAM(5)/1/
     61  DATA IPARAM(6),IPARAM(7),IPARAM(8),IPARAM(9)/0,0,0,0/
     62  INTEGER :: J4SAVE,IWHICH,IVALUE
     63  !***FIRST EXECUTABLE STATEMENT  J4SAVE
     64  J4SAVE = IPARAM(IWHICH)
     65  IF (ISET) IPARAM(IWHICH) = IVALUE
     66  RETURN
     67END FUNCTION J4SAVE

  • LMDZ6/trunk/libf/misc/juldate.f90

    r5245 r5246  
    22! $Id$
    33!
    4         subroutine juldate(ian,imoi,ijou,oh,om,os,tjd,tjdsec)
    5 c       Sous-routine de changement de date:
    6 c       gregorien>>>date julienne
    7 c       En entree:an,mois,jour,heure,min.,sec.
    8 c       En sortie:tjd
    9         IMPLICIT NONE
    10         INTEGER,INTENT(IN) :: ian,imoi,ijou,oh,om,os
    11         REAL,INTENT(OUT) :: tjd,tjdsec
    12        
    13         REAL frac,year,rmon,cf,a,b
    14         INTEGER ojou
    15        
    16         frac=((os/60.+om)/60.+oh)/24.
    17         ojou=dble(ijou)+frac
    18             year=dble(ian)
    19             rmon=dble(imoi)
    20         if (imoi .le. 2) then
    21             year=year-1.
    22             rmon=rmon+12.
    23         endif
    24         cf=year+(rmon/100.)+(ojou/10000.)
    25         if (cf .ge. 1582.1015) then
    26             a=int(year/100)
    27             b=2-a+int(a/4)
    28         else
    29             b=0
    30         endif
    31         tjd=int(365.25*year)+int(30.6001*(rmon+1))+int(ojou)
    32      +   +1720994.5+b
    33         tjdsec=(ojou-int(ojou))+(tjd-int(tjd))
    34         tjd=int(tjd)+int(tjdsec)
    35         tjdsec=tjdsec-int(tjdsec)
    36         return
    37         end
     4  subroutine juldate(ian,imoi,ijou,oh,om,os,tjd,tjdsec)
     5  ! Sous-routine de changement de date:
     6  ! gregorien>>>date julienne
     7  ! En entree:an,mois,jour,heure,min.,sec.
     8  ! En sortie:tjd
     9    IMPLICIT NONE
     10    INTEGER,INTENT(IN) :: ian,imoi,ijou,oh,om,os
     11    REAL,INTENT(OUT) :: tjd,tjdsec
     12
     13    REAL :: frac,year,rmon,cf,a,b
     14    INTEGER :: ojou
     15
     16    frac=((os/60.+om)/60.+oh)/24.
     17    ojou=dble(ijou)+frac
     18        year=dble(ian)
     19        rmon=dble(imoi)
     20    if (imoi .le. 2) then
     21        year=year-1.
     22        rmon=rmon+12.
     23    endif
     24    cf=year+(rmon/100.)+(ojou/10000.)
     25    if (cf .ge. 1582.1015) then
     26        a=int(year/100)
     27        b=2-a+int(a/4)
     28    else
     29        b=0
     30    endif
     31    tjd=int(365.25*year)+int(30.6001*(rmon+1))+int(ojou) &
     32          +1720994.5+b
     33    tjdsec=(ojou-int(ojou))+(tjd-int(tjd))
     34    tjd=int(tjd)+int(tjdsec)
     35    tjdsec=tjdsec-int(tjdsec)
     36    return
     37end subroutine juldate
    3838
    3939

  • LMDZ6/trunk/libf/misc/minmax.f90

    r5245 r5246  
    22! $Header$
    33!
    4        SUBROUTINE minmax(imax, xi, zmin, zmax )
    5 c
    6 c      P. Le Van
     4 SUBROUTINE minmax(imax, xi, zmin, zmax )
     5  !
     6  !  P. Le Van
    77
    8        INTEGER imax
    9        REAL    xi(imax)
    10        REAL    zmin,zmax
    11        INTEGER i
     8   INTEGER :: imax
     9   REAL :: xi(imax)
     10   REAL :: zmin,zmax
     11   INTEGER :: i
    1212
    13        zmin = xi(1)
    14        zmax = xi(1)
     13   zmin = xi(1)
     14   zmax = xi(1)
    1515
    16        DO i = 2, imax
    17          zmin = MIN( zmin,xi(i) )
    18          zmax = MAX( zmax,xi(i) )
    19        ENDDO
     16   DO i = 2, imax
     17     zmin = MIN( zmin,xi(i) )
     18     zmax = MAX( zmax,xi(i) )
     19   ENDDO
    2020
    21        RETURN
    22        END
     21   RETURN
     22END SUBROUTINE minmax
    2323

  • LMDZ6/trunk/libf/misc/minmax2.f90

    r5245 r5246  
    22! $Header$
    33!
    4        SUBROUTINE minmax2(imax, jmax, lmax, xi, zmin, zmax )
    5 c
    6        INTEGER lmax,jmax,imax
    7        REAL xi(imax*jmax*lmax)
    8        REAL zmin,zmax
    9        INTEGER i
    10    
    11        zmin = xi(1)
    12        zmax = xi(1)
     4 SUBROUTINE minmax2(imax, jmax, lmax, xi, zmin, zmax )
     5  !
     6   INTEGER :: lmax,jmax,imax
     7   REAL :: xi(imax*jmax*lmax)
     8   REAL :: zmin,zmax
     9   INTEGER :: i
    1310
    14        DO i = 2, imax*jmax*lmax
    15          zmin = MIN( zmin,xi(i) )
    16          zmax = MAX( zmax,xi(i) )
    17        ENDDO
     11   zmin = xi(1)
     12   zmax = xi(1)
    1813
    19        RETURN
    20        END
     14   DO i = 2, imax*jmax*lmax
     15     zmin = MIN( zmin,xi(i) )
     16     zmax = MAX( zmax,xi(i) )
     17   ENDDO
     18
     19   RETURN
     20END SUBROUTINE minmax2

  • LMDZ6/trunk/libf/misc/pchdf.f90

    r5245 r5246  
    1 *DECK PCHDF
    2       REAL FUNCTION PCHDF (K, X, S, IERR)
    3 C***BEGIN PROLOGUE  PCHDF
    4 C***SUBSIDIARY
    5 C***PURPOSE  Computes divided differences for PCHCE and PCHSP
    6 C***LIBRARY   SLATEC (PCHIP)
    7 C***TYPE      SINGLE PRECISION (PCHDF-S, DPCHDF-D)
    8 C***AUTHOR  Fritsch, F. N., (LLNL)
    9 C***DESCRIPTION
    10 C
    11 C          PCHDF:   PCHIP Finite Difference Formula
    12 C
    13 C    Uses a divided difference formulation to compute a K-point approx-
    14 C    imation to the derivative at X(K) based on the data in X and S.
    15 C
    16 C    Called by  PCHCE  and  PCHSP  to compute 3- and 4-point boundary
    17 C    derivative approximations.
    18 C
    19 C ----------------------------------------------------------------------
    20 C
    21 C    On input:
    22 C        K      is the order of the desired derivative approximation.
    23 C               K must be at least 3 (error return if not).
    24 C        X      contains the K values of the independent variable.
    25 C               X need not be ordered, but the values **MUST** be
    26 C               distinct.  (Not checked here.)
    27 C        S      contains the associated slope values:
    28 C                  S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1.
    29 C               (Note that S need only be of length K-1.)
    30 C
    31 C    On return:
    32 C        S      will be destroyed.
    33 C        IERR   will be set to -1 if K.LT.2 .
    34 C        PCHDF  will be set to the desired derivative approximation if
    35 C               IERR=0 or to zero if IERR=-1.
    36 C
    37 C ----------------------------------------------------------------------
    38 C
    39 C***SEE ALSO  PCHCE, PCHSP
    40 C***REFERENCES  Carl de Boor, A Practical Guide to Splines, Springer-
    41 C                 Verlag, New York, 1978, pp. 10-16.
    42 C***ROUTINES CALLED  XERMSG
    43 C***REVISION HISTORY  (YYMMDD)
    44 C   820503  DATE WRITTEN
    45 C   820805  Converted to SLATEC library version.
    46 C   870813  Minor cosmetic changes.
    47 C   890411  Added SAVE statements (Vers. 3.2).
    48 C   890411  REVISION DATE from Version 3.2
    49 C   891214  Prologue converted to Version 4.0 format.  (BAB)
    50 C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
    51 C   900328  Added TYPE section.  (WRB)
    52 C   910408  Updated AUTHOR and DATE WRITTEN sections in prologue.  (WRB)
    53 C   920429  Revised format and order of references.  (WRB,FNF)
    54 C   930503  Improved purpose.  (FNF)
    55 C***END PROLOGUE  PCHDF
    56 C
    57 C**End
    58 C
    59 C  DECLARE ARGUMENTS.
    60 C
    61       INTEGER K, IERR
    62       REAL X(K), S(K)
    63 C
    64 C  DECLARE LOCAL VARIABLES.
    65 C
    66       INTEGER I, J
    67       REAL VALUE, ZERO
    68       SAVE ZERO
    69       DATA  ZERO /0./
    70 C
    71 C  CHECK FOR LEGAL VALUE OF K.
    72 C
    73 C***FIRST EXECUTABLE STATEMENT  PCHDF
    74       IF (K .LT. 3)  GO TO 5001
    75 C
    76 C  COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL.
    77 C
    78       DO 10  J = 2, K-1
    79          DO 9  I = 1, K-J
    80             S(I) = (S(I+1)-S(I))/(X(I+J)-X(I))
    81     9    CONTINUE
    82    10 CONTINUE
    83 C
    84 C  EVALUATE DERIVATIVE AT X(K).
    85 C
    86       VALUE = S(1)
    87       DO 20  I = 2, K-1
    88          VALUE = S(I) + VALUE*(X(K)-X(I))
    89    20 CONTINUE
    90 C
    91 C  NORMAL RETURN.
    92 C
    93       IERR = 0
    94       PCHDF = VALUE
    95       RETURN
    96 C
    97 C  ERROR RETURN.
    98 C
    99  5001 CONTINUE
    100 C    K.LT.3 RETURN.
    101       IERR = -1
    102       CALL XERMSG ('SLATEC', 'PCHDF', 'K LESS THAN THREE', IERR, 1)
    103       PCHDF = ZERO
    104       RETURN
    105 C------------- LAST LINE OF PCHDF FOLLOWS ------------------------------
    106       END
     1!DECK PCHDF
     2REAL FUNCTION PCHDF (K, X, S, IERR)
     3  !***BEGIN PROLOGUE  PCHDF
     4  !***SUBSIDIARY
     5  !***PURPOSE  Computes divided differences for PCHCE and PCHSP
     6  !***LIBRARY   SLATEC (PCHIP)
     7  !***TYPE      SINGLE PRECISION (PCHDF-S, DPCHDF-D)
     8  !***AUTHOR  Fritsch, F. N., (LLNL)
     9  !***DESCRIPTION
     10  !
     11  !      PCHDF:   PCHIP Finite Difference Formula
     12  !
     13  ! Uses a divided difference formulation to compute a K-point approx-
     14  ! imation to the derivative at X(K) based on the data in X and S.
     15  !
     16  ! Called by  PCHCE  and  PCHSP  to compute 3- and 4-point boundary
     17  ! derivative approximations.
     18  !
     19  ! ----------------------------------------------------------------------
     20  !
     21  ! On input:
     22  !    K      is the order of the desired derivative approximation.
     23  !           K must be at least 3 (error return if not).
     24  !    X      contains the K values of the independent variable.
     25  !           X need not be ordered, but the values **MUST** be
     26  !           distinct.  (Not checked here.)
     27  !    S      contains the associated slope values:
     28  !              S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1.
     29  !           (Note that S need only be of length K-1.)
     30  !
     31  ! On return:
     32  !    S      will be destroyed.
     33  !    IERR   will be set to -1 if K.LT.2 .
     34  !    PCHDF  will be set to the desired derivative approximation if
     35  !           IERR=0 or to zero if IERR=-1.
     36  !
     37  ! ----------------------------------------------------------------------
     38  !
     39  !***SEE ALSO  PCHCE, PCHSP
     40  !***REFERENCES  Carl de Boor, A Practical Guide to Splines, Springer-
     41  !             Verlag, New York, 1978, pp. 10-16.
     42  !***ROUTINES CALLED  XERMSG
     43  !***REVISION HISTORY  (YYMMDD)
     44  !   820503  DATE WRITTEN
     45  !   820805  Converted to SLATEC library version.
     46  !   870813  Minor cosmetic changes.
     47  !   890411  Added SAVE statements (Vers. 3.2).
     48  !   890411  REVISION DATE from Version 3.2
     49  !   891214  Prologue converted to Version 4.0 format.  (BAB)
     50  !   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
     51  !   900328  Added TYPE section.  (WRB)
     52  !   910408  Updated AUTHOR and DATE WRITTEN sections in prologue.  (WRB)
     53  !   920429  Revised format and order of references.  (WRB,FNF)
     54  !   930503  Improved purpose.  (FNF)
     55  !***END PROLOGUE  PCHDF
     56  !
     57  !**End
     58  !
     59  !  DECLARE ARGUMENTS.
     60  !
     61  INTEGER :: K, IERR
     62  REAL :: X(K), S(K)
     63  !
     64  !  DECLARE LOCAL VARIABLES.
     65  !
     66  INTEGER :: I, J
     67  REAL :: VALUE, ZERO
     68  SAVE ZERO
     69  DATA  ZERO /0./
     70  !
     71  !  CHECK FOR LEGAL VALUE OF K.
     72  !
     73  !***FIRST EXECUTABLE STATEMENT  PCHDF
     74  IF (K .LT. 3)  GO TO 5001
     75  !
     76  !  COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL.
     77  !
     78  DO  J = 2, K-1
     79     DO  I = 1, K-J
     80        S(I) = (S(I+1)-S(I))/(X(I+J)-X(I))
     81     END DO
     82  END DO
     83  !
     84  !  EVALUATE DERIVATIVE AT X(K).
     85  !
     86  VALUE = S(1)
     87  DO  I = 2, K-1
     88     VALUE = S(I) + VALUE*(X(K)-X(I))
     89  END DO
     90  !
     91  !  NORMAL RETURN.
     92  !
     93  IERR = 0
     94  PCHDF = VALUE
     95  RETURN
     96  !
     97  !  ERROR RETURN.
     98  !
     99 5001   CONTINUE
     100  ! K.LT.3 RETURN.
     101  IERR = -1
     102  CALL XERMSG ('SLATEC', 'PCHDF', 'K LESS THAN THREE', IERR, 1)
     103  PCHDF = ZERO
     104  RETURN
     105  !------------- LAST LINE OF PCHDF FOLLOWS ------------------------------
     106END FUNCTION PCHDF

  • LMDZ6/trunk/libf/misc/pchfe.f90

    r5245 r5246  
    1 *DECK PCHFE
    2       SUBROUTINE PCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)
    3 C***BEGIN PROLOGUE  PCHFE
    4 C***PURPOSE  Evaluate a piecewise cubic Hermite function at an array of
    5 C            points.  May be used by itself for Hermite interpolation,
    6 C            or as an evaluator for PCHIM or PCHIC.
    7 C***LIBRARY   SLATEC (PCHIP)
    8 C***CATEGORY  E3
    9 C***TYPE      SINGLE PRECISION (PCHFE-S, DPCHFE-D)
    10 C***KEYWORDS  CUBIC HERMITE EVALUATION, HERMITE INTERPOLATION, PCHIP,
    11 C             PIECEWISE CUBIC EVALUATION
    12 C***AUTHOR  Fritsch, F. N., (LLNL)
    13 C             Lawrence Livermore National Laboratory
    14 C             P.O. Box 808  (L-316)
    15 C             Livermore, CA  94550
    16 C             FTS 532-4275, (510) 422-4275
    17 C***DESCRIPTION
    18 C
    19 C          PCHFE:  Piecewise Cubic Hermite Function Evaluator
    20 C
    21 C    Evaluates the cubic Hermite function defined by  N, X, F, D  at
    22 C    the points  XE(J), J=1(1)NE.
    23 C
    24 C    To provide compatibility with PCHIM and PCHIC, includes an
    25 C    increment between successive values of the F- and D-arrays.
    26 C
    27 C ----------------------------------------------------------------------
    28 C
    29 C  Calling sequence:
    30 C
    31 C        PARAMETER  (INCFD = ...)
    32 C        INTEGER  N, NE, IERR
    33 C        REAL  X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE)
    34 C        LOGICAL  SKIP
    35 C
    36 C        CALL  PCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)
    37 C
    38 C   Parameters:
    39 C
    40 C    N -- (input) number of data points.  (Error return if N.LT.2 .)
    41 C
    42 C    X -- (input) real array of independent variable values.  The
    43 C           elements of X must be strictly increasing:
    44 C                X(I-1) .LT. X(I),  I = 2(1)N.
    45 C           (Error return if not.)
    46 C
    47 C    F -- (input) real array of function values.  F(1+(I-1)*INCFD) is
    48 C           the value corresponding to X(I).
    49 C
    50 C    D -- (input) real array of derivative values.  D(1+(I-1)*INCFD) is
    51 C           the value corresponding to X(I).
    52 C
    53 C    INCFD -- (input) increment between successive values in F and D.
    54 C           (Error return if  INCFD.LT.1 .)
    55 C
    56 C    SKIP -- (input/output) logical variable which should be set to
    57 C           .TRUE. if the user wishes to skip checks for validity of
    58 C           preceding parameters, or to .FALSE. otherwise.
    59 C           This will save time in case these checks have already
    60 C           been performed (say, in PCHIM or PCHIC).
    61 C           SKIP will be set to .TRUE. on normal return.
    62 C
    63 C    NE -- (input) number of evaluation points.  (Error return if
    64 C           NE.LT.1 .)
    65 C
    66 C    XE -- (input) real array of points at which the function is to be
    67 C           evaluated.
    68 C
    69 C          NOTES:
    70 C           1. The evaluation will be most efficient if the elements
    71 C              of XE are increasing relative to X;
    72 C              that is,   XE(J) .GE. X(I)
    73 C              implies    XE(K) .GE. X(I),  all K.GE.J .
    74 C           2. If any of the XE are outside the interval [X(1),X(N)],
    75 C              values are extrapolated from the nearest extreme cubic,
    76 C              and a warning error is returned.
    77 C
    78 C    FE -- (output) real array of values of the cubic Hermite function
    79 C           defined by  N, X, F, D  at the points  XE.
    80 C
    81 C    IERR -- (output) error flag.
    82 C           Normal return:
    83 C              IERR = 0  (no errors).
    84 C           Warning error:
    85 C              IERR.GT.0  means that extrapolation was performed at
    86 C                 IERR points.
    87 C           "Recoverable" errors:
    88 C              IERR = -1  if N.LT.2 .
    89 C              IERR = -2  if INCFD.LT.1 .
    90 C              IERR = -3  if the X-array is not strictly increasing.
    91 C              IERR = -4  if NE.LT.1 .
    92 C             (The FE-array has not been changed in any of these cases.)
    93 C               NOTE:  The above errors are checked in the order listed,
    94 C                   and following arguments have **NOT** been validated.
    95 C
    96 C***REFERENCES  (NONE)
    97 C***ROUTINES CALLED  CHFEV, XERMSG
    98 C***REVISION HISTORY  (YYMMDD)
    99 C   811020  DATE WRITTEN
    100 C   820803  Minor cosmetic changes for release 1.
    101 C   870707  Minor cosmetic changes to prologue.
    102 C   890531  Changed all specific intrinsics to generic.  (WRB)
    103 C   890831  Modified array declarations.  (WRB)
    104 C   890831  REVISION DATE from Version 3.2
    105 C   891214  Prologue converted to Version 4.0 format.  (BAB)
    106 C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
    107 C***END PROLOGUE  PCHFE
    108 C  Programming notes:
    109 C
    110 C    1. To produce a double precision version, simply:
    111 C        a. Change PCHFE to DPCHFE, and CHFEV to DCHFEV, wherever they
    112 C           occur,
    113 C        b. Change the real declaration to double precision,
    114 C
    115 C    2. Most of the coding between the call to CHFEV and the end of
    116 C        the IR-loop could be eliminated if it were permissible to
    117 C        assume that XE is ordered relative to X.
    118 C
    119 C    3. CHFEV does not assume that X1 is less than X2.  thus, it would
    120 C        be possible to write a version of PCHFE that assumes a strict-
    121 C        ly decreasing X-array by simply running the IR-loop backwards
    122 C        (and reversing the order of appropriate tests).
    123 C
    124 C    4. The present code has a minor bug, which I have decided is not
    125 C        worth the effort that would be required to fix it.
    126 C        If XE contains points in [X(N-1),X(N)], followed by points .LT.
    127 C        X(N-1), followed by points .GT.X(N), the extrapolation points
    128 C        will be counted (at least) twice in the total returned in IERR.
    129 C
    130 C  DECLARE ARGUMENTS.
    131 C
    132       INTEGER N, INCFD, NE, IERR
    133       REAL X(*), F(INCFD,*), D(INCFD,*), XE(*), FE(*)
    134       LOGICAL SKIP
    135 C
    136 C  DECLARE LOCAL VARIABLES.
    137 C
    138       INTEGER I, IERC, IR, J, JFIRST, NEXT(2), NJ
    139 C
    140 C  VALIDITY-CHECK ARGUMENTS.
    141 C
    142 C***FIRST EXECUTABLE STATEMENT  PCHFE
    143       IF (SKIP)  GO TO 5
    144 C
    145       IF ( N.LT.2 )  GO TO 5001
    146       IF ( INCFD.LT.1 )  GO TO 5002
    147       DO 1  I = 2, N
    148          IF ( X(I).LE.X(I-1) )  GO TO 5003
    149     1 CONTINUE
    150 C
    151 C  FUNCTION DEFINITION IS OK, GO ON.
    152 C
    153     5 CONTINUE
    154       IF ( NE.LT.1 )  GO TO 5004
    155       IERR = 0
    156       SKIP = .TRUE.
    157 C
    158 C  LOOP OVER INTERVALS.        (   INTERVAL INDEX IS  IL = IR-1  . )
    159 C                              ( INTERVAL IS X(IL).LE.X.LT.X(IR) . )
    160       JFIRST = 1
    161       IR = 2
    162    10 CONTINUE
    163 C
    164 C    SKIP OUT OF LOOP IF HAVE PROCESSED ALL EVALUATION POINTS.
    165 C
    166          IF (JFIRST .GT. NE)  GO TO 5000
    167 C
    168 C    LOCATE ALL POINTS IN INTERVAL.
    169 C
    170          DO 20  J = JFIRST, NE
    171             IF (XE(J) .GE. X(IR))  GO TO 30
    172    20    CONTINUE
    173          J = NE + 1
    174          GO TO 40
    175 C
    176 C    HAVE LOCATED FIRST POINT BEYOND INTERVAL.
    177 C
    178    30    CONTINUE
    179          IF (IR .EQ. N)  J = NE + 1
    180 C
    181    40    CONTINUE
    182          NJ = J - JFIRST
    183 C
    184 C    SKIP EVALUATION IF NO POINTS IN INTERVAL.
    185 C
    186          IF (NJ .EQ. 0)  GO TO 50
    187 C
    188 C    EVALUATE CUBIC AT XE(I),  I = JFIRST (1) J-1 .
    189 C
    190 C       ----------------------------------------------------------------
    191         CALL CHFEV (X(IR-1),X(IR), F(1,IR-1),F(1,IR), D(1,IR-1),D(1,IR),
    192      *              NJ, XE(JFIRST), FE(JFIRST), NEXT, IERC)
    193 C      ----------------------------------------------------------------
    194          IF (IERC .LT. 0)  GO TO 5005
    195 C
    196          IF (NEXT(2) .EQ. 0)  GO TO 42
    197 C        IF (NEXT(2) .GT. 0)  THEN
    198 C           IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(2) TO THE
    199 C           RIGHT OF X(IR).
    200 C
    201             IF (IR .LT. N)  GO TO 41
    202 C          IF (IR .EQ. N)  THEN
    203 C              THESE ARE ACTUALLY EXTRAPOLATION POINTS.
    204                IERR = IERR + NEXT(2)
    205                GO TO 42
    206    41       CONTINUE
    207 C          ELSE
    208 C              WE SHOULD NEVER HAVE GOTTEN HERE.
    209                GO TO 5005
    210 C          ENDIF
    211 C        ENDIF
    212    42    CONTINUE
    213 C
    214          IF (NEXT(1) .EQ. 0)  GO TO 49
    215 C        IF (NEXT(1) .GT. 0)  THEN
    216 C           IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(1) TO THE
    217 C           LEFT OF X(IR-1).
    218 C
    219             IF (IR .GT. 2)  GO TO 43
    220 C          IF (IR .EQ. 2)  THEN
    221 C              THESE ARE ACTUALLY EXTRAPOLATION POINTS.
    222                IERR = IERR + NEXT(1)
    223                GO TO 49
    224    43       CONTINUE
    225 C          ELSE
    226 C              XE IS NOT ORDERED RELATIVE TO X, SO MUST ADJUST
    227 C              EVALUATION INTERVAL.
    228 C
    229 C              FIRST, LOCATE FIRST POINT TO LEFT OF X(IR-1).
    230                DO 44  I = JFIRST, J-1
    231                   IF (XE(I) .LT. X(IR-1))  GO TO 45
    232    44          CONTINUE
    233 C              NOTE-- CANNOT DROP THROUGH HERE UNLESS THERE IS AN ERROR
    234 C                     IN CHFEV.
    235                GO TO 5005
    236 C
    237    45          CONTINUE
    238 C              RESET J.  (THIS WILL BE THE NEW JFIRST.)
    239                J = I
    240 C
    241 C              NOW FIND OUT HOW FAR TO BACK UP IN THE X-ARRAY.
    242                DO 46  I = 1, IR-1
    243                   IF (XE(J) .LT. X(I)) GO TO 47
    244    46          CONTINUE
    245 C              NB-- CAN NEVER DROP THROUGH HERE, SINCE XE(J).LT.X(IR-1).
    246 C
    247    47          CONTINUE
    248 C              AT THIS POINT, EITHER  XE(J) .LT. X(1)
    249 C                 OR      X(I-1) .LE. XE(J) .LT. X(I) .
    250 C              RESET IR, RECOGNIZING THAT IT WILL BE INCREMENTED BEFORE
    251 C              CYCLING.
    252                IR = MAX(1, I-1)
    253 C          ENDIF
    254 C        ENDIF
    255    49    CONTINUE
    256 C
    257          JFIRST = J
    258 C
    259 C    END OF IR-LOOP.
    260 C
    261    50 CONTINUE
    262       IR = IR + 1
    263       IF (IR .LE. N)  GO TO 10
    264 C
    265 C  NORMAL RETURN.
    266 C
    267  5000 CONTINUE
    268       RETURN
    269 C
    270 C  ERROR RETURNS.
    271 C
    272  5001 CONTINUE
    273 C    N.LT.2 RETURN.
    274       IERR = -1
    275       CALL XERMSG ('SLATEC', 'PCHFE',
    276      +   'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
    277       RETURN
    278 C
    279  5002 CONTINUE
    280 C    INCFD.LT.1 RETURN.
    281       IERR = -2
    282       CALL XERMSG ('SLATEC', 'PCHFE', 'INCREMENT LESS THAN ONE', IERR,
    283      +   1)
    284       RETURN
    285 C
    286  5003 CONTINUE
    287 C    X-ARRAY NOT STRICTLY INCREASING.
    288       IERR = -3
    289       CALL XERMSG ('SLATEC', 'PCHFE', 'X-ARRAY NOT STRICTLY INCREASING'
    290      +   , IERR, 1)
    291       RETURN
    292 C
    293  5004 CONTINUE
    294 C    NE.LT.1 RETURN.
    295       IERR = -4
    296       CALL XERMSG ('SLATEC', 'PCHFE',
    297      +   'NUMBER OF EVALUATION POINTS LESS THAN ONE', IERR, 1)
    298       RETURN
    299 C
    300  5005 CONTINUE
    301 C    ERROR RETURN FROM CHFEV.
    302 C   *** THIS CASE SHOULD NEVER OCCUR ***
    303       IERR = -5
    304       CALL XERMSG ('SLATEC', 'PCHFE',
    305      +   'ERROR RETURN FROM CHFEV -- FATAL', IERR, 2)
    306       RETURN
    307 C------------- LAST LINE OF PCHFE FOLLOWS ------------------------------
    308       END
     1!DECK PCHFE
     2SUBROUTINE PCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)
     3  !***BEGIN PROLOGUE  PCHFE
     4  !***PURPOSE  Evaluate a piecewise cubic Hermite function at an array of
     5         ! points.  May be used by itself for Hermite interpolation,
     6         ! or as an evaluator for PCHIM or PCHIC.
     7  !***LIBRARY   SLATEC (PCHIP)
     8  !***CATEGORY  E3
     9  !***TYPE      SINGLE PRECISION (PCHFE-S, DPCHFE-D)
     10  !***KEYWORDS  CUBIC HERMITE EVALUATION, HERMITE INTERPOLATION, PCHIP,
     11         !  PIECEWISE CUBIC EVALUATION
     12  !***AUTHOR  Fritsch, F. N., (LLNL)
     13         !  Lawrence Livermore National Laboratory
     14         !  P.O. Box 808  (L-316)
     15         !  Livermore, CA  94550
     16         !  FTS 532-4275, (510) 422-4275
     17  !***DESCRIPTION
     18  !
     19  !      PCHFE:  Piecewise Cubic Hermite Function Evaluator
     20  !
     21  ! Evaluates the cubic Hermite function defined by  N, X, F, D  at
     22  ! the points  XE(J), J=1(1)NE.
     23  !
     24  ! To provide compatibility with PCHIM and PCHIC, includes an
     25  ! increment between successive values of the F- and D-arrays.
     26  !
     27  ! ----------------------------------------------------------------------
     28  !
     29  !  Calling sequence:
     30  !
     31  !    PARAMETER  (INCFD = ...)
     32  !    INTEGER  N, NE, IERR
     33  !    REAL  X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE)
     34  !    LOGICAL  SKIP
     35  !
     36  !    CALL  PCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)
     37  !
     38  !   Parameters:
     39  !
     40  ! N -- (input) number of data points.  (Error return if N.LT.2 .)
     41  !
     42  ! X -- (input) real array of independent variable values.  The
     43  !       elements of X must be strictly increasing:
     44  !            X(I-1) .LT. X(I),  I = 2(1)N.
     45  !       (Error return if not.)
     46  !
     47  ! F -- (input) real array of function values.  F(1+(I-1)*INCFD) is
     48  !       the value corresponding to X(I).
     49  !
     50  ! D -- (input) real array of derivative values.  D(1+(I-1)*INCFD) is
     51  !       the value corresponding to X(I).
     52  !
     53  ! INCFD -- (input) increment between successive values in F and D.
     54  !       (Error return if  INCFD.LT.1 .)
     55  !
     56  ! SKIP -- (input/output) logical variable which should be set to
     57  !       .TRUE. if the user wishes to skip checks for validity of
     58  !       preceding parameters, or to .FALSE. otherwise.
     59  !       This will save time in case these checks have already
     60  !       been performed (say, in PCHIM or PCHIC).
     61  !       SKIP will be set to .TRUE. on normal return.
     62  !
     63  ! NE -- (input) number of evaluation points.  (Error return if
     64  !       NE.LT.1 .)
     65  !
     66  ! XE -- (input) real array of points at which the function is to be
     67  !       evaluated.
     68  !
     69  !      NOTES:
     70  !       1. The evaluation will be most efficient if the elements
     71  !          of XE are increasing relative to X;
     72  !          that is,   XE(J) .GE. X(I)
     73  !          implies    XE(K) .GE. X(I),  all K.GE.J .
     74  !       2. If any of the XE are outside the interval [X(1),X(N)],
     75  !          values are extrapolated from the nearest extreme cubic,
     76  !          and a warning error is returned.
     77  !
     78  ! FE -- (output) real array of values of the cubic Hermite function
     79  !       defined by  N, X, F, D  at the points  XE.
     80  !
     81  ! IERR -- (output) error flag.
     82  !       Normal return:
     83  !          IERR = 0  (no errors).
     84  !       Warning error:
     85  !          IERR.GT.0  means that extrapolation was performed at
     86  !             IERR points.
     87  !       "Recoverable" errors:
     88  !          IERR = -1  if N.LT.2 .
     89  !          IERR = -2  if INCFD.LT.1 .
     90  !          IERR = -3  if the X-array is not strictly increasing.
     91  !          IERR = -4  if NE.LT.1 .
     92  !         (The FE-array has not been changed in any of these cases.)
     93  !           NOTE:  The above errors are checked in the order listed,
     94  !               and following arguments have **NOT** been validated.
     95  !
     96  !***REFERENCES  (NONE)
     97  !***ROUTINES CALLED  CHFEV, XERMSG
     98  !***REVISION HISTORY  (YYMMDD)
     99  !   811020  DATE WRITTEN
     100  !   820803  Minor cosmetic changes for release 1.
     101  !   870707  Minor cosmetic changes to prologue.
     102  !   890531  Changed all specific intrinsics to generic.  (WRB)
     103  !   890831  Modified array declarations.  (WRB)
     104  !   890831  REVISION DATE from Version 3.2
     105  !   891214  Prologue converted to Version 4.0 format.  (BAB)
     106  !   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
     107  !***END PROLOGUE  PCHFE
     108  !  Programming notes:
     109  !
     110  ! 1. To produce a double precision version, simply:
     111  !    a. Change PCHFE to DPCHFE, and CHFEV to DCHFEV, wherever they
     112  !       occur,
     113  !    b. Change the real declaration to double precision,
     114  !
     115  ! 2. Most of the coding between the call to CHFEV and the end of
     116  !    the IR-loop could be eliminated if it were permissible to
     117  !    assume that XE is ordered relative to X.
     118  !
     119  ! 3. CHFEV does not assume that X1 is less than X2.  thus, it would
     120  !    be possible to write a version of PCHFE that assumes a strict-
     121  !    ly decreasing X-array by simply running the IR-loop backwards
     122  !    (and reversing the order of appropriate tests).
     123  !
     124  ! 4. The present code has a minor bug, which I have decided is not
     125  !    worth the effort that would be required to fix it.
     126  !    If XE contains points in [X(N-1),X(N)], followed by points .LT.
     127  !    X(N-1), followed by points .GT.X(N), the extrapolation points
     128  !    will be counted (at least) twice in the total returned in IERR.
     129  !
     130  !  DECLARE ARGUMENTS.
     131  !
     132  INTEGER :: N, INCFD, NE, IERR
     133  REAL :: X(*), F(INCFD,*), D(INCFD,*), XE(*), FE(*)
     134  LOGICAL :: SKIP
     135  !
     136  !  DECLARE LOCAL VARIABLES.
     137  !
     138  INTEGER :: I, IERC, IR, J, JFIRST, NEXT(2), NJ
     139  !
     140  !  VALIDITY-CHECK ARGUMENTS.
     141  !
     142  !***FIRST EXECUTABLE STATEMENT  PCHFE
     143  IF (SKIP)  GO TO 5
     144  !
     145  IF ( N.LT.2 )  GO TO 5001
     146  IF ( INCFD.LT.1 )  GO TO 5002
     147  DO  I = 2, N
     148     IF ( X(I).LE.X(I-1) )  GO TO 5003
     149  END DO
     150  !
     151  !  FUNCTION DEFINITION IS OK, GO ON.
     152  !
     153    5   CONTINUE
     154  IF ( NE.LT.1 )  GO TO 5004
     155  IERR = 0
     156  SKIP = .TRUE.
     157  !
     158  !  LOOP OVER INTERVALS.        (   INTERVAL INDEX IS  IL = IR-1  . )
     159  !                              ( INTERVAL IS X(IL).LE.X.LT.X(IR) . )
     160  JFIRST = 1
     161  IR = 2
     162   10   CONTINUE
     163  !
     164  ! SKIP OUT OF LOOP IF HAVE PROCESSED ALL EVALUATION POINTS.
     165  !
     166     IF (JFIRST .GT. NE)  GO TO 5000
     167  !
     168  ! LOCATE ALL POINTS IN INTERVAL.
     169  !
     170     DO  J = JFIRST, NE
     171        IF (XE(J) .GE. X(IR))  GO TO 30
     172     END DO
     173     J = NE + 1
     174     GO TO 40
     175  !
     176  ! HAVE LOCATED FIRST POINT BEYOND INTERVAL.
     177  !
     178   30   CONTINUE
     179     IF (IR .EQ. N)  J = NE + 1
     180  !
     181   40   CONTINUE
     182     NJ = J - JFIRST
     183  !
     184  ! SKIP EVALUATION IF NO POINTS IN INTERVAL.
     185  !
     186     IF (NJ .EQ. 0)  GO TO 50
     187  !
     188  ! EVALUATE CUBIC AT XE(I),  I = JFIRST (1) J-1 .
     189  !
     190  !   ----------------------------------------------------------------
     191    CALL CHFEV (X(IR-1),X(IR), F(1,IR-1),F(1,IR), D(1,IR-1),D(1,IR), &
     192          NJ, XE(JFIRST), FE(JFIRST), NEXT, IERC)
     193    ! ----------------------------------------------------------------
     194     IF (IERC .LT. 0)  GO TO 5005
     195  !
     196     IF (NEXT(2) .EQ. 0)  GO TO 42
     197     ! IF (NEXT(2) .GT. 0)  THEN
     198     !    IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(2) TO THE
     199     !    RIGHT OF X(IR).
     200  !
     201        IF (IR .LT. N)  GO TO 41
     202        ! IF (IR .EQ. N)  THEN
     203        !    THESE ARE ACTUALLY EXTRAPOLATION POINTS.
     204           IERR = IERR + NEXT(2)
     205           GO TO 42
     206   41    CONTINUE
     207        ! ELSE
     208        !    WE SHOULD NEVER HAVE GOTTEN HERE.
     209           GO TO 5005
     210        ! ENDIF
     211     ! ENDIF
     212   42   CONTINUE
     213  !
     214     IF (NEXT(1) .EQ. 0)  GO TO 49
     215     ! IF (NEXT(1) .GT. 0)  THEN
     216     !    IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(1) TO THE
     217     !    LEFT OF X(IR-1).
     218  !
     219        IF (IR .GT. 2)  GO TO 43
     220        ! IF (IR .EQ. 2)  THEN
     221        !    THESE ARE ACTUALLY EXTRAPOLATION POINTS.
     222           IERR = IERR + NEXT(1)
     223           GO TO 49
     224   43    CONTINUE
     225        ! ELSE
     226        !    XE IS NOT ORDERED RELATIVE TO X, SO MUST ADJUST
     227        !    EVALUATION INTERVAL.
     228  !
     229  !          FIRST, LOCATE FIRST POINT TO LEFT OF X(IR-1).
     230           DO  I = JFIRST, J-1
     231              IF (XE(I) .LT. X(IR-1))  GO TO 45
     232           END DO
     233           ! NOTE-- CANNOT DROP THROUGH HERE UNLESS THERE IS AN ERROR
     234           !        IN CHFEV.
     235           GO TO 5005
     236  !
     237   45       CONTINUE
     238           ! RESET J.  (THIS WILL BE THE NEW JFIRST.)
     239           J = I
     240  !
     241  !          NOW FIND OUT HOW FAR TO BACK UP IN THE X-ARRAY.
     242           DO  I = 1, IR-1
     243              IF (XE(J) .LT. X(I)) GO TO 47
     244           END DO
     245           ! NB-- CAN NEVER DROP THROUGH HERE, SINCE XE(J).LT.X(IR-1).
     246  !
     247   47       CONTINUE
     248           ! AT THIS POINT, EITHER  XE(J) .LT. X(1)
     249           !    OR      X(I-1) .LE. XE(J) .LT. X(I) .
     250           ! RESET IR, RECOGNIZING THAT IT WILL BE INCREMENTED BEFORE
     251           ! CYCLING.
     252           IR = MAX(1, I-1)
     253        ! ENDIF
     254     ! ENDIF
     255   49   CONTINUE
     256  !
     257     JFIRST = J
     258  !
     259  ! END OF IR-LOOP.
     260  !
     261   50   CONTINUE
     262  IR = IR + 1
     263  IF (IR .LE. N)  GO TO 10
     264  !
     265  !  NORMAL RETURN.
     266  !
     267 5000   CONTINUE
     268  RETURN
     269  !
     270  !  ERROR RETURNS.
     271  !
     272 5001   CONTINUE
     273  ! N.LT.2 RETURN.
     274  IERR = -1
     275  CALL XERMSG ('SLATEC', 'PCHFE', &
     276        'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
     277  RETURN
     278  !
     279 5002   CONTINUE
     280  ! INCFD.LT.1 RETURN.
     281  IERR = -2
     282  CALL XERMSG ('SLATEC', 'PCHFE', 'INCREMENT LESS THAN ONE', IERR, &
     283        1)
     284  RETURN
     285  !
     286 5003   CONTINUE
     287  ! X-ARRAY NOT STRICTLY INCREASING.
     288  IERR = -3
     289  CALL XERMSG ('SLATEC', 'PCHFE', 'X-ARRAY NOT STRICTLY INCREASING' &
     290        , IERR, 1)
     291  RETURN
     292  !
     293 5004   CONTINUE
     294  ! NE.LT.1 RETURN.
     295  IERR = -4
     296  CALL XERMSG ('SLATEC', 'PCHFE', &
     297        'NUMBER OF EVALUATION POINTS LESS THAN ONE', IERR, 1)
     298  RETURN
     299  !
     300 5005   CONTINUE
     301  ! ERROR RETURN FROM CHFEV.
     302  !   *** THIS CASE SHOULD NEVER OCCUR ***
     303  IERR = -5
     304  CALL XERMSG ('SLATEC', 'PCHFE', &
     305        'ERROR RETURN FROM CHFEV -- FATAL', IERR, 2)
     306  RETURN
     307  !------------- LAST LINE OF PCHFE FOLLOWS ------------------------------
     308END SUBROUTINE PCHFE

  • LMDZ6/trunk/libf/misc/pchsp.f90

    r5245 r5246  
    1 *DECK PCHSP
    2       SUBROUTINE PCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR)
    3 C***BEGIN PROLOGUE  PCHSP
    4 C***PURPOSE  Set derivatives needed to determine the Hermite represen-
    5 C            tation of the cubic spline interpolant to given data, with
    6 C            specified boundary conditions.
    7 C***LIBRARY   SLATEC (PCHIP)
    8 C***CATEGORY  E1A
    9 C***TYPE      SINGLE PRECISION (PCHSP-S, DPCHSP-D)
    10 C***KEYWORDS  CUBIC HERMITE INTERPOLATION, PCHIP,
    11 C             PIECEWISE CUBIC INTERPOLATION, SPLINE INTERPOLATION
    12 C***AUTHOR  Fritsch, F. N., (LLNL)
    13 C             Lawrence Livermore National Laboratory
    14 C             P.O. Box 808  (L-316)
    15 C             Livermore, CA  94550
    16 C             FTS 532-4275, (510) 422-4275
    17 C***DESCRIPTION
    18 C
    19 C          PCHSP:   Piecewise Cubic Hermite Spline
    20 C
    21 C    Computes the Hermite representation of the cubic spline inter-
    22 C    polant to the data given in X and F satisfying the boundary
    23 C    conditions specified by IC and VC.
    24 C
    25 C    To facilitate two-dimensional applications, includes an increment
    26 C    between successive values of the F- and D-arrays.
    27 C
    28 C    The resulting piecewise cubic Hermite function may be evaluated
    29 C    by PCHFE or PCHFD.
    30 C
    31 C    NOTE:  This is a modified version of C. de Boor's cubic spline
    32 C            routine CUBSPL.
    33 C
    34 C ----------------------------------------------------------------------
    35 C
    36 C  Calling sequence:
    37 C
    38 C        PARAMETER  (INCFD = ...)
    39 C        INTEGER  IC(2), N, NWK, IERR
    40 C        REAL  VC(2), X(N), F(INCFD,N), D(INCFD,N), WK(NWK)
    41 C
    42 C        CALL  PCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR)
    43 C
    44 C   Parameters:
    45 C
    46 C    IC -- (input) integer array of length 2 specifying desired
    47 C           boundary conditions:
    48 C           IC(1) = IBEG, desired condition at beginning of data.
    49 C           IC(2) = IEND, desired condition at end of data.
    50 C
    51 C           IBEG = 0  to set D(1) so that the third derivative is con-
    52 C              tinuous at X(2).  This is the "not a knot" condition
    53 C              provided by de Boor's cubic spline routine CUBSPL.
    54 C              < This is the default boundary condition. >
    55 C           IBEG = 1  if first derivative at X(1) is given in VC(1).
    56 C           IBEG = 2  if second derivative at X(1) is given in VC(1).
    57 C           IBEG = 3  to use the 3-point difference formula for D(1).
    58 C                     (Reverts to the default b.c. if N.LT.3 .)
    59 C           IBEG = 4  to use the 4-point difference formula for D(1).
    60 C                     (Reverts to the default b.c. if N.LT.4 .)
    61 C          NOTES:
    62 C           1. An error return is taken if IBEG is out of range.
    63 C           2. For the "natural" boundary condition, use IBEG=2 and
    64 C              VC(1)=0.
    65 C
    66 C           IEND may take on the same values as IBEG, but applied to
    67 C           derivative at X(N).  In case IEND = 1 or 2, the value is
    68 C           given in VC(2).
    69 C
    70 C          NOTES:
    71 C           1. An error return is taken if IEND is out of range.
    72 C           2. For the "natural" boundary condition, use IEND=2 and
    73 C              VC(2)=0.
    74 C
    75 C    VC -- (input) real array of length 2 specifying desired boundary
    76 C           values, as indicated above.
    77 C           VC(1) need be set only if IC(1) = 1 or 2 .
    78 C           VC(2) need be set only if IC(2) = 1 or 2 .
    79 C
    80 C    N -- (input) number of data points.  (Error return if N.LT.2 .)
    81 C
    82 C    X -- (input) real array of independent variable values.  The
    83 C           elements of X must be strictly increasing:
    84 C                X(I-1) .LT. X(I),  I = 2(1)N.
    85 C           (Error return if not.)
    86 C
    87 C    F -- (input) real array of dependent variable values to be inter-
    88 C           polated.  F(1+(I-1)*INCFD) is value corresponding to X(I).
    89 C
    90 C    D -- (output) real array of derivative values at the data points.
    91 C           These values will determine the cubic spline interpolant
    92 C           with the requested boundary conditions.
    93 C           The value corresponding to X(I) is stored in
    94 C                D(1+(I-1)*INCFD),  I=1(1)N.
    95 C           No other entries in D are changed.
    96 C
    97 C    INCFD -- (input) increment between successive values in F and D.
    98 C           This argument is provided primarily for 2-D applications.
    99 C           (Error return if  INCFD.LT.1 .)
    100 C
    101 C    WK -- (scratch) real array of working storage.
    102 C
    103 C    NWK -- (input) length of work array.
    104 C           (Error return if NWK.LT.2*N .)
    105 C
    106 C    IERR -- (output) error flag.
    107 C           Normal return:
    108 C              IERR = 0  (no errors).
    109 C           "Recoverable" errors:
    110 C              IERR = -1  if N.LT.2 .
    111 C              IERR = -2  if INCFD.LT.1 .
    112 C              IERR = -3  if the X-array is not strictly increasing.
    113 C              IERR = -4  if IBEG.LT.0 or IBEG.GT.4 .
    114 C              IERR = -5  if IEND.LT.0 of IEND.GT.4 .
    115 C              IERR = -6  if both of the above are true.
    116 C              IERR = -7  if NWK is too small.
    117 C               NOTE:  The above errors are checked in the order listed,
    118 C                   and following arguments have **NOT** been validated.
    119 C             (The D-array has not been changed in any of these cases.)
    120 C              IERR = -8  in case of trouble solving the linear system
    121 C                         for the interior derivative values.
    122 C             (The D-array may have been changed in this case.)
    123 C             (             Do **NOT** use it!                )
    124 C
    125 C***REFERENCES  Carl de Boor, A Practical Guide to Splines, Springer-
    126 C                 Verlag, New York, 1978, pp. 53-59.
    127 C***ROUTINES CALLED  PCHDF, XERMSG
    128 C***REVISION HISTORY  (YYMMDD)
    129 C   820503  DATE WRITTEN
    130 C   820804  Converted to SLATEC library version.
    131 C   870707  Minor cosmetic changes to prologue.
    132 C   890411  Added SAVE statements (Vers. 3.2).
    133 C   890703  Corrected category record.  (WRB)
    134 C   890831  Modified array declarations.  (WRB)
    135 C   890831  REVISION DATE from Version 3.2
    136 C   891214  Prologue converted to Version 4.0 format.  (BAB)
    137 C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
    138 C   920429  Revised format and order of references.  (WRB,FNF)
    139 C***END PROLOGUE  PCHSP
    140 C  Programming notes:
    141 C
    142 C    To produce a double precision version, simply:
    143 C        a. Change PCHSP to DPCHSP wherever it occurs,
    144 C        b. Change the real declarations to double precision, and
    145 C        c. Change the constants ZERO, HALF, ... to double precision.
    146 C
    147 C  DECLARE ARGUMENTS.
    148 C
    149       INTEGER IC(2), N, INCFD, NWK, IERR
    150       REAL VC(2), X(*), F(INCFD,*), D(INCFD,*), WK(2,*)
    151 C
    152 C  DECLARE LOCAL VARIABLES.
    153 C
    154       INTEGER IBEG, IEND, INDEX, J, NM1
    155       REAL G, HALF, ONE, STEMP(3), THREE, TWO, XTEMP(4), ZERO
    156       SAVE ZERO, HALF, ONE, TWO, THREE
    157       REAL PCHDF
    158 C
    159       DATA  ZERO /0./,  HALF /0.5/,  ONE /1./,  TWO /2./,  THREE /3./
    160 C
    161 C  VALIDITY-CHECK ARGUMENTS.
    162 C
    163 C***FIRST EXECUTABLE STATEMENT  PCHSP
    164       IF ( N.LT.2 )  GO TO 5001
    165       IF ( INCFD.LT.1 )  GO TO 5002
    166       DO 1  J = 2, N
    167          IF ( X(J).LE.X(J-1) )  GO TO 5003
    168     1 CONTINUE
    169 C
    170       IBEG = IC(1)
    171       IEND = IC(2)
    172       IERR = 0
    173       IF ( (IBEG.LT.0).OR.(IBEG.GT.4) )  IERR = IERR - 1
    174       IF ( (IEND.LT.0).OR.(IEND.GT.4) )  IERR = IERR - 2
    175       IF ( IERR.LT.0 )  GO TO 5004
    176 C
    177 C  FUNCTION DEFINITION IS OK -- GO ON.
    178 C
    179       IF ( NWK .LT. 2*N )  GO TO 5007
    180 C
    181 C  COMPUTE FIRST DIFFERENCES OF X SEQUENCE AND STORE IN WK(1,.). ALSO,
    182 C  COMPUTE FIRST DIVIDED DIFFERENCE OF DATA AND STORE IN WK(2,.).
    183       DO 5  J=2,N
    184          WK(1,J) = X(J) - X(J-1)
    185          WK(2,J) = (F(1,J) - F(1,J-1))/WK(1,J)
    186     5 CONTINUE
    187 C
    188 C  SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL.
    189 C
    190       IF ( IBEG.GT.N )  IBEG = 0
    191       IF ( IEND.GT.N )  IEND = 0
    192 C
    193 C  SET UP FOR BOUNDARY CONDITIONS.
    194 C
    195       IF ( (IBEG.EQ.1).OR.(IBEG.EQ.2) )  THEN
    196          D(1,1) = VC(1)
    197       ELSE IF (IBEG .GT. 2)  THEN
    198 C        PICK UP FIRST IBEG POINTS, IN REVERSE ORDER.
    199          DO 10  J = 1, IBEG
    200             INDEX = IBEG-J+1
    201 C          INDEX RUNS FROM IBEG DOWN TO 1.
    202             XTEMP(J) = X(INDEX)
    203             IF (J .LT. IBEG)  STEMP(J) = WK(2,INDEX)
    204    10    CONTINUE
    205 C                --------------------------------
    206          D(1,1) = PCHDF (IBEG, XTEMP, STEMP, IERR)
    207 C                --------------------------------
    208          IF (IERR .NE. 0)  GO TO 5009
    209          IBEG = 1
    210       ENDIF
    211 C
    212       IF ( (IEND.EQ.1).OR.(IEND.EQ.2) )  THEN
    213          D(1,N) = VC(2)
    214       ELSE IF (IEND .GT. 2)  THEN
    215 C        PICK UP LAST IEND POINTS.
    216          DO 15  J = 1, IEND
    217             INDEX = N-IEND+J
    218 C          INDEX RUNS FROM N+1-IEND UP TO N.
    219             XTEMP(J) = X(INDEX)
    220             IF (J .LT. IEND)  STEMP(J) = WK(2,INDEX+1)
    221    15    CONTINUE
    222 C                --------------------------------
    223          D(1,N) = PCHDF (IEND, XTEMP, STEMP, IERR)
    224 C                --------------------------------
    225          IF (IERR .NE. 0)  GO TO 5009
    226          IEND = 1
    227       ENDIF
    228 C
    229 C --------------------( BEGIN CODING FROM CUBSPL )--------------------
    230 C
    231 C  **** A TRIDIAGONAL LINEAR SYSTEM FOR THE UNKNOWN SLOPES S(J) OF
    232 C  F  AT X(J), J=1,...,N, IS GENERATED AND THEN SOLVED BY GAUSS ELIM-
    233 C  INATION, WITH S(J) ENDING UP IN D(1,J), ALL J.
    234 C    WK(1,.) AND WK(2,.) ARE USED FOR TEMPORARY STORAGE.
    235 C
    236 C  CONSTRUCT FIRST EQUATION FROM FIRST BOUNDARY CONDITION, OF THE FORM
    237 C             WK(2,1)*S(1) + WK(1,1)*S(2) = D(1,1)
    238 C
    239       IF (IBEG .EQ. 0)  THEN
    240          IF (N .EQ. 2)  THEN
    241 C          NO CONDITION AT LEFT END AND N = 2.
    242             WK(2,1) = ONE
    243             WK(1,1) = ONE
    244             D(1,1) = TWO*WK(2,2)
    245          ELSE
    246 C          NOT-A-KNOT CONDITION AT LEFT END AND N .GT. 2.
    247             WK(2,1) = WK(1,3)
    248             WK(1,1) = WK(1,2) + WK(1,3)
    249             D(1,1) =((WK(1,2) + TWO*WK(1,1))*WK(2,2)*WK(1,3)
    250      *                        + WK(1,2)**2*WK(2,3)) / WK(1,1)
    251          ENDIF
    252       ELSE IF (IBEG .EQ. 1)  THEN
    253 C        SLOPE PRESCRIBED AT LEFT END.
    254          WK(2,1) = ONE
    255          WK(1,1) = ZERO
    256       ELSE
    257 C        SECOND DERIVATIVE PRESCRIBED AT LEFT END.
    258          WK(2,1) = TWO
    259          WK(1,1) = ONE
    260          D(1,1) = THREE*WK(2,2) - HALF*WK(1,2)*D(1,1)
    261       ENDIF
    262 C
    263 C  IF THERE ARE INTERIOR KNOTS, GENERATE THE CORRESPONDING EQUATIONS AND
    264 C  CARRY OUT THE FORWARD PASS OF GAUSS ELIMINATION, AFTER WHICH THE J-TH
    265 C  EQUATION READS    WK(2,J)*S(J) + WK(1,J)*S(J+1) = D(1,J).
    266 C
    267       NM1 = N-1
    268       IF (NM1 .GT. 1)  THEN
    269          DO 20 J=2,NM1
    270             IF (WK(2,J-1) .EQ. ZERO)  GO TO 5008
    271             G = -WK(1,J+1)/WK(2,J-1)
    272             D(1,J) = G*D(1,J-1)
    273      *                  + THREE*(WK(1,J)*WK(2,J+1) + WK(1,J+1)*WK(2,J))
    274             WK(2,J) = G*WK(1,J-1) + TWO*(WK(1,J) + WK(1,J+1))
    275    20    CONTINUE
    276       ENDIF
    277 C
    278 C  CONSTRUCT LAST EQUATION FROM SECOND BOUNDARY CONDITION, OF THE FORM
    279 C           (-G*WK(2,N-1))*S(N-1) + WK(2,N)*S(N) = D(1,N)
    280 C
    281 C    IF SLOPE IS PRESCRIBED AT RIGHT END, ONE CAN GO DIRECTLY TO BACK-
    282 C    SUBSTITUTION, SINCE ARRAYS HAPPEN TO BE SET UP JUST RIGHT FOR IT
    283 C    AT THIS POINT.
    284       IF (IEND .EQ. 1)  GO TO 30
    285 C
    286       IF (IEND .EQ. 0)  THEN
    287          IF (N.EQ.2 .AND. IBEG.EQ.0)  THEN
    288 C          NOT-A-KNOT AT RIGHT ENDPOINT AND AT LEFT ENDPOINT AND N = 2.
    289             D(1,2) = WK(2,2)
    290             GO TO 30
    291          ELSE IF ((N.EQ.2) .OR. (N.EQ.3 .AND. IBEG.EQ.0))  THEN
    292 C          EITHER (N=3 AND NOT-A-KNOT ALSO AT LEFT) OR (N=2 AND *NOT*
    293 C          NOT-A-KNOT AT LEFT END POINT).
    294             D(1,N) = TWO*WK(2,N)
    295             WK(2,N) = ONE
    296             IF (WK(2,N-1) .EQ. ZERO)  GO TO 5008
    297             G = -ONE/WK(2,N-1)
    298          ELSE
    299 C          NOT-A-KNOT AND N .GE. 3, AND EITHER N.GT.3 OR  ALSO NOT-A-
    300 C          KNOT AT LEFT END POINT.
    301             G = WK(1,N-1) + WK(1,N)
    302 C          DO NOT NEED TO CHECK FOLLOWING DENOMINATORS (X-DIFFERENCES).
    303             D(1,N) = ((WK(1,N)+TWO*G)*WK(2,N)*WK(1,N-1)
    304      *                  + WK(1,N)**2*(F(1,N-1)-F(1,N-2))/WK(1,N-1))/G
    305             IF (WK(2,N-1) .EQ. ZERO)  GO TO 5008
    306             G = -G/WK(2,N-1)
    307             WK(2,N) = WK(1,N-1)
    308          ENDIF
    309       ELSE
    310 C        SECOND DERIVATIVE PRESCRIBED AT RIGHT ENDPOINT.
    311          D(1,N) = THREE*WK(2,N) + HALF*WK(1,N)*D(1,N)
    312          WK(2,N) = TWO
    313          IF (WK(2,N-1) .EQ. ZERO)  GO TO 5008
    314          G = -ONE/WK(2,N-1)
    315       ENDIF
    316 C
    317 C  COMPLETE FORWARD PASS OF GAUSS ELIMINATION.
    318 C
    319       WK(2,N) = G*WK(1,N-1) + WK(2,N)
    320       IF (WK(2,N) .EQ. ZERO)   GO TO 5008
    321       D(1,N) = (G*D(1,N-1) + D(1,N))/WK(2,N)
    322 C
    323 C  CARRY OUT BACK SUBSTITUTION
    324 C
    325    30 CONTINUE
    326       DO 40 J=NM1,1,-1
    327          IF (WK(2,J) .EQ. ZERO)  GO TO 5008
    328          D(1,J) = (D(1,J) - WK(1,J)*D(1,J+1))/WK(2,J)
    329    40 CONTINUE
    330 C --------------------(  END  CODING FROM CUBSPL )--------------------
    331 C
    332 C  NORMAL RETURN.
    333 C
    334       RETURN
    335 C
    336 C  ERROR RETURNS.
    337 C
    338  5001 CONTINUE
    339 C    N.LT.2 RETURN.
    340       IERR = -1
    341       CALL XERMSG ('SLATEC', 'PCHSP',
    342      +   'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
    343       RETURN
    344 C
    345  5002 CONTINUE
    346 C    INCFD.LT.1 RETURN.
    347       IERR = -2
    348       CALL XERMSG ('SLATEC', 'PCHSP', 'INCREMENT LESS THAN ONE', IERR,
    349      +   1)
    350       RETURN
    351 C
    352  5003 CONTINUE
    353 C    X-ARRAY NOT STRICTLY INCREASING.
    354       IERR = -3
    355       CALL XERMSG ('SLATEC', 'PCHSP', 'X-ARRAY NOT STRICTLY INCREASING'
    356      +   , IERR, 1)
    357       RETURN
    358 C
    359  5004 CONTINUE
    360 C    IC OUT OF RANGE RETURN.
    361       IERR = IERR - 3
    362       CALL XERMSG ('SLATEC', 'PCHSP', 'IC OUT OF RANGE', IERR, 1)
    363       RETURN
    364 C
    365  5007 CONTINUE
    366 C    NWK TOO SMALL RETURN.
    367       IERR = -7
    368       CALL XERMSG ('SLATEC', 'PCHSP', 'WORK ARRAY TOO SMALL', IERR, 1)
    369       RETURN
    370 C
    371  5008 CONTINUE
    372 C    SINGULAR SYSTEM.
    373 C   *** THEORETICALLY, THIS CAN ONLY OCCUR IF SUCCESSIVE X-VALUES   ***
    374 C   *** ARE EQUAL, WHICH SHOULD ALREADY HAVE BEEN CAUGHT (IERR=-3). ***
    375       IERR = -8
    376       CALL XERMSG ('SLATEC', 'PCHSP', 'SINGULAR LINEAR SYSTEM', IERR,
    377      +   1)
    378       RETURN
    379 C
    380  5009 CONTINUE
    381 C    ERROR RETURN FROM PCHDF.
    382 C   *** THIS CASE SHOULD NEVER OCCUR ***
    383       IERR = -9
    384       CALL XERMSG ('SLATEC', 'PCHSP', 'ERROR RETURN FROM PCHDF', IERR,
    385      +   1)
    386       RETURN
    387 C------------- LAST LINE OF PCHSP FOLLOWS ------------------------------
    388       END
     1!DECK PCHSP
     2SUBROUTINE PCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR)
     3  !***BEGIN PROLOGUE  PCHSP
     4  !***PURPOSE  Set derivatives needed to determine the Hermite represen-
     5         ! tation of the cubic spline interpolant to given data, with
     6         ! specified boundary conditions.
     7  !***LIBRARY   SLATEC (PCHIP)
     8  !***CATEGORY  E1A
     9  !***TYPE      SINGLE PRECISION (PCHSP-S, DPCHSP-D)
     10  !***KEYWORDS  CUBIC HERMITE INTERPOLATION, PCHIP,
     11         !  PIECEWISE CUBIC INTERPOLATION, SPLINE INTERPOLATION
     12  !***AUTHOR  Fritsch, F. N., (LLNL)
     13         !  Lawrence Livermore National Laboratory
     14         !  P.O. Box 808  (L-316)
     15         !  Livermore, CA  94550
     16         !  FTS 532-4275, (510) 422-4275
     17  !***DESCRIPTION
     18  !
     19  !      PCHSP:   Piecewise Cubic Hermite Spline
     20  !
     21  ! Computes the Hermite representation of the cubic spline inter-
     22  ! polant to the data given in X and F satisfying the boundary
     23  ! conditions specified by IC and VC.
     24  !
     25  ! To facilitate two-dimensional applications, includes an increment
     26  ! between successive values of the F- and D-arrays.
     27  !
     28  ! The resulting piecewise cubic Hermite function may be evaluated
     29  ! by PCHFE or PCHFD.
     30  !
     31  ! NOTE:  This is a modified version of C. de Boor's cubic spline
     32  !        routine CUBSPL.
     33  !
     34  ! ----------------------------------------------------------------------
     35  !
     36  !  Calling sequence:
     37  !
     38  !    PARAMETER  (INCFD = ...)
     39  !    INTEGER  IC(2), N, NWK, IERR
     40  !    REAL  VC(2), X(N), F(INCFD,N), D(INCFD,N), WK(NWK)
     41  !
     42  !    CALL  PCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR)
     43  !
     44  !   Parameters:
     45  !
     46  ! IC -- (input) integer array of length 2 specifying desired
     47  !       boundary conditions:
     48  !       IC(1) = IBEG, desired condition at beginning of data.
     49  !       IC(2) = IEND, desired condition at end of data.
     50  !
     51  !       IBEG = 0  to set D(1) so that the third derivative is con-
     52  !          tinuous at X(2).  This is the "not a knot" condition
     53  !          provided by de Boor's cubic spline routine CUBSPL.
     54  !          < This is the default boundary condition. >
     55  !       IBEG = 1  if first derivative at X(1) is given in VC(1).
     56  !       IBEG = 2  if second derivative at X(1) is given in VC(1).
     57  !       IBEG = 3  to use the 3-point difference formula for D(1).
     58  !                 (Reverts to the default b.c. if N.LT.3 .)
     59  !       IBEG = 4  to use the 4-point difference formula for D(1).
     60  !                 (Reverts to the default b.c. if N.LT.4 .)
     61  !      NOTES:
     62  !       1. An error return is taken if IBEG is out of range.
     63  !       2. For the "natural" boundary condition, use IBEG=2 and
     64  !          VC(1)=0.
     65  !
     66  !       IEND may take on the same values as IBEG, but applied to
     67  !       derivative at X(N).  In case IEND = 1 or 2, the value is
     68  !       given in VC(2).
     69  !
     70  !      NOTES:
     71  !       1. An error return is taken if IEND is out of range.
     72  !       2. For the "natural" boundary condition, use IEND=2 and
     73  !          VC(2)=0.
     74  !
     75  ! VC -- (input) real array of length 2 specifying desired boundary
     76  !       values, as indicated above.
     77  !       VC(1) need be set only if IC(1) = 1 or 2 .
     78  !       VC(2) need be set only if IC(2) = 1 or 2 .
     79  !
     80  ! N -- (input) number of data points.  (Error return if N.LT.2 .)
     81  !
     82  ! X -- (input) real array of independent variable values.  The
     83  !       elements of X must be strictly increasing:
     84  !            X(I-1) .LT. X(I),  I = 2(1)N.
     85  !       (Error return if not.)
     86  !
     87  ! F -- (input) real array of dependent variable values to be inter-
     88  !       polated.  F(1+(I-1)*INCFD) is value corresponding to X(I).
     89  !
     90  ! D -- (output) real array of derivative values at the data points.
     91  !       These values will determine the cubic spline interpolant
     92  !       with the requested boundary conditions.
     93  !       The value corresponding to X(I) is stored in
     94  !            D(1+(I-1)*INCFD),  I=1(1)N.
     95  !       No other entries in D are changed.
     96  !
     97  ! INCFD -- (input) increment between successive values in F and D.
     98  !       This argument is provided primarily for 2-D applications.
     99  !       (Error return if  INCFD.LT.1 .)
     100  !
     101  ! WK -- (scratch) real array of working storage.
     102  !
     103  ! NWK -- (input) length of work array.
     104  !       (Error return if NWK.LT.2*N .)
     105  !
     106  ! IERR -- (output) error flag.
     107  !       Normal return:
     108  !          IERR = 0  (no errors).
     109  !       "Recoverable" errors:
     110  !          IERR = -1  if N.LT.2 .
     111  !          IERR = -2  if INCFD.LT.1 .
     112  !          IERR = -3  if the X-array is not strictly increasing.
     113  !          IERR = -4  if IBEG.LT.0 or IBEG.GT.4 .
     114  !          IERR = -5  if IEND.LT.0 of IEND.GT.4 .
     115  !          IERR = -6  if both of the above are true.
     116  !          IERR = -7  if NWK is too small.
     117  !           NOTE:  The above errors are checked in the order listed,
     118  !               and following arguments have **NOT** been validated.
     119  !         (The D-array has not been changed in any of these cases.)
     120  !          IERR = -8  in case of trouble solving the linear system
     121  !                     for the interior derivative values.
     122  !         (The D-array may have been changed in this case.)
     123  !         (             Do **NOT** use it!                )
     124  !
     125  !***REFERENCES  Carl de Boor, A Practical Guide to Splines, Springer-
     126  !             Verlag, New York, 1978, pp. 53-59.
     127  !***ROUTINES CALLED  PCHDF, XERMSG
     128  !***REVISION HISTORY  (YYMMDD)
     129  !   820503  DATE WRITTEN
     130  !   820804  Converted to SLATEC library version.
     131  !   870707  Minor cosmetic changes to prologue.
     132  !   890411  Added SAVE statements (Vers. 3.2).
     133  !   890703  Corrected category record.  (WRB)
     134  !   890831  Modified array declarations.  (WRB)
     135  !   890831  REVISION DATE from Version 3.2
     136  !   891214  Prologue converted to Version 4.0 format.  (BAB)
     137  !   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
     138  !   920429  Revised format and order of references.  (WRB,FNF)
     139  !***END PROLOGUE  PCHSP
     140  !  Programming notes:
     141  !
     142  ! To produce a double precision version, simply:
     143  !    a. Change PCHSP to DPCHSP wherever it occurs,
     144  !    b. Change the real declarations to double precision, and
     145  !    c. Change the constants ZERO, HALF, ... to double precision.
     146  !
     147  !  DECLARE ARGUMENTS.
     148  !
     149  INTEGER :: IC(2), N, INCFD, NWK, IERR
     150  REAL :: VC(2), X(*), F(INCFD,*), D(INCFD,*), WK(2,*)
     151  !
     152  !  DECLARE LOCAL VARIABLES.
     153  !
     154  INTEGER :: IBEG, IEND, INDEX, J, NM1
     155  REAL :: G, HALF, ONE, STEMP(3), THREE, TWO, XTEMP(4), ZERO
     156  SAVE ZERO, HALF, ONE, TWO, THREE
     157  REAL :: PCHDF
     158  !
     159  DATA  ZERO /0./,  HALF /0.5/,  ONE /1./,  TWO /2./,  THREE /3./
     160  !
     161  !  VALIDITY-CHECK ARGUMENTS.
     162  !
     163  !***FIRST EXECUTABLE STATEMENT  PCHSP
     164  IF ( N.LT.2 )  GO TO 5001
     165  IF ( INCFD.LT.1 )  GO TO 5002
     166  DO  J = 2, N
     167     IF ( X(J).LE.X(J-1) )  GO TO 5003
     168  END DO
     169  !
     170  IBEG = IC(1)
     171  IEND = IC(2)
     172  IERR = 0
     173  IF ( (IBEG.LT.0).OR.(IBEG.GT.4) )  IERR = IERR - 1
     174  IF ( (IEND.LT.0).OR.(IEND.GT.4) )  IERR = IERR - 2
     175  IF ( IERR.LT.0 )  GO TO 5004
     176  !
     177  !  FUNCTION DEFINITION IS OK -- GO ON.
     178  !
     179  IF ( NWK .LT. 2*N )  GO TO 5007
     180  !
     181  !  COMPUTE FIRST DIFFERENCES OF X SEQUENCE AND STORE IN WK(1,.). ALSO,
     182  !  COMPUTE FIRST DIVIDED DIFFERENCE OF DATA AND STORE IN WK(2,.).
     183  DO  J=2,N
     184     WK(1,J) = X(J) - X(J-1)
     185     WK(2,J) = (F(1,J) - F(1,J-1))/WK(1,J)
     186  END DO
     187  !
     188  !  SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL.
     189  !
     190  IF ( IBEG.GT.N )  IBEG = 0
     191  IF ( IEND.GT.N )  IEND = 0
     192  !
     193  !  SET UP FOR BOUNDARY CONDITIONS.
     194  !
     195  IF ( (IBEG.EQ.1).OR.(IBEG.EQ.2) )  THEN
     196     D(1,1) = VC(1)
     197  ELSE IF (IBEG .GT. 2)  THEN
     198     ! PICK UP FIRST IBEG POINTS, IN REVERSE ORDER.
     199     DO  J = 1, IBEG
     200        INDEX = IBEG-J+1
     201        ! INDEX RUNS FROM IBEG DOWN TO 1.
     202        XTEMP(J) = X(INDEX)
     203        IF (J .LT. IBEG)  STEMP(J) = WK(2,INDEX)
     204     END DO
     205              ! --------------------------------
     206     D(1,1) = PCHDF (IBEG, XTEMP, STEMP, IERR)
     207              ! --------------------------------
     208     IF (IERR .NE. 0)  GO TO 5009
     209     IBEG = 1
     210  ENDIF
     211  !
     212  IF ( (IEND.EQ.1).OR.(IEND.EQ.2) )  THEN
     213     D(1,N) = VC(2)
     214  ELSE IF (IEND .GT. 2)  THEN
     215     ! PICK UP LAST IEND POINTS.
     216     DO  J = 1, IEND
     217        INDEX = N-IEND+J
     218        ! INDEX RUNS FROM N+1-IEND UP TO N.
     219        XTEMP(J) = X(INDEX)
     220        IF (J .LT. IEND)  STEMP(J) = WK(2,INDEX+1)
     221     END DO
     222              ! --------------------------------
     223     D(1,N) = PCHDF (IEND, XTEMP, STEMP, IERR)
     224              ! --------------------------------
     225     IF (IERR .NE. 0)  GO TO 5009
     226     IEND = 1
     227  ENDIF
     228  !
     229  ! --------------------( BEGIN CODING FROM CUBSPL )--------------------
     230  !
     231  !  **** A TRIDIAGONAL LINEAR SYSTEM FOR THE UNKNOWN SLOPES S(J) OF
     232  !  F  AT X(J), J=1,...,N, IS GENERATED AND THEN SOLVED BY GAUSS ELIM-
     233  !  INATION, WITH S(J) ENDING UP IN D(1,J), ALL J.
     234  ! WK(1,.) AND WK(2,.) ARE USED FOR TEMPORARY STORAGE.
     235  !
     236  !  CONSTRUCT FIRST EQUATION FROM FIRST BOUNDARY CONDITION, OF THE FORM
     237  !         WK(2,1)*S(1) + WK(1,1)*S(2) = D(1,1)
     238  !
     239  IF (IBEG .EQ. 0)  THEN
     240     IF (N .EQ. 2)  THEN
     241        ! NO CONDITION AT LEFT END AND N = 2.
     242        WK(2,1) = ONE
     243        WK(1,1) = ONE
     244        D(1,1) = TWO*WK(2,2)
     245     ELSE
     246        ! NOT-A-KNOT CONDITION AT LEFT END AND N .GT. 2.
     247        WK(2,1) = WK(1,3)
     248        WK(1,1) = WK(1,2) + WK(1,3)
     249        D(1,1) =((WK(1,2) + TWO*WK(1,1))*WK(2,2)*WK(1,3) &
     250              + WK(1,2)**2*WK(2,3)) / WK(1,1)
     251     ENDIF
     252  ELSE IF (IBEG .EQ. 1)  THEN
     253     ! SLOPE PRESCRIBED AT LEFT END.
     254     WK(2,1) = ONE
     255     WK(1,1) = ZERO
     256  ELSE
     257     ! SECOND DERIVATIVE PRESCRIBED AT LEFT END.
     258     WK(2,1) = TWO
     259     WK(1,1) = ONE
     260     D(1,1) = THREE*WK(2,2) - HALF*WK(1,2)*D(1,1)
     261  ENDIF
     262  !
     263  !  IF THERE ARE INTERIOR KNOTS, GENERATE THE CORRESPONDING EQUATIONS AND
     264  !  CARRY OUT THE FORWARD PASS OF GAUSS ELIMINATION, AFTER WHICH THE J-TH
     265  !  EQUATION READS    WK(2,J)*S(J) + WK(1,J)*S(J+1) = D(1,J).
     266  !
     267  NM1 = N-1
     268  IF (NM1 .GT. 1)  THEN
     269     DO J=2,NM1
     270        IF (WK(2,J-1) .EQ. ZERO)  GO TO 5008
     271        G = -WK(1,J+1)/WK(2,J-1)
     272        D(1,J) = G*D(1,J-1) &
     273              + THREE*(WK(1,J)*WK(2,J+1) + WK(1,J+1)*WK(2,J))
     274        WK(2,J) = G*WK(1,J-1) + TWO*(WK(1,J) + WK(1,J+1))
     275     END DO
     276  ENDIF
     277  !
     278  !  CONSTRUCT LAST EQUATION FROM SECOND BOUNDARY CONDITION, OF THE FORM
     279  !       (-G*WK(2,N-1))*S(N-1) + WK(2,N)*S(N) = D(1,N)
     280  !
     281  ! IF SLOPE IS PRESCRIBED AT RIGHT END, ONE CAN GO DIRECTLY TO BACK-
     282  ! SUBSTITUTION, SINCE ARRAYS HAPPEN TO BE SET UP JUST RIGHT FOR IT
     283  ! AT THIS POINT.
     284  IF (IEND .EQ. 1)  GO TO 30
     285  !
     286  IF (IEND .EQ. 0)  THEN
     287     IF (N.EQ.2 .AND. IBEG.EQ.0)  THEN
     288        ! NOT-A-KNOT AT RIGHT ENDPOINT AND AT LEFT ENDPOINT AND N = 2.
     289        D(1,2) = WK(2,2)
     290        GO TO 30
     291     ELSE IF ((N.EQ.2) .OR. (N.EQ.3 .AND. IBEG.EQ.0))  THEN
     292        ! EITHER (N=3 AND NOT-A-KNOT ALSO AT LEFT) OR (N=2 AND *NOT*
     293        ! NOT-A-KNOT AT LEFT END POINT).
     294        D(1,N) = TWO*WK(2,N)
     295        WK(2,N) = ONE
     296        IF (WK(2,N-1) .EQ. ZERO)  GO TO 5008
     297        G = -ONE/WK(2,N-1)
     298     ELSE
     299        ! NOT-A-KNOT AND N .GE. 3, AND EITHER N.GT.3 OR  ALSO NOT-A-
     300        ! KNOT AT LEFT END POINT.
     301        G = WK(1,N-1) + WK(1,N)
     302        ! DO NOT NEED TO CHECK FOLLOWING DENOMINATORS (X-DIFFERENCES).
     303        D(1,N) = ((WK(1,N)+TWO*G)*WK(2,N)*WK(1,N-1) &
     304              + WK(1,N)**2*(F(1,N-1)-F(1,N-2))/WK(1,N-1))/G
     305        IF (WK(2,N-1) .EQ. ZERO)  GO TO 5008
     306        G = -G/WK(2,N-1)
     307        WK(2,N) = WK(1,N-1)
     308     ENDIF
     309  ELSE
     310     ! SECOND DERIVATIVE PRESCRIBED AT RIGHT ENDPOINT.
     311     D(1,N) = THREE*WK(2,N) + HALF*WK(1,N)*D(1,N)
     312     WK(2,N) = TWO
     313     IF (WK(2,N-1) .EQ. ZERO)  GO TO 5008
     314     G = -ONE/WK(2,N-1)
     315  ENDIF
     316  !
     317  !  COMPLETE FORWARD PASS OF GAUSS ELIMINATION.
     318  !
     319  WK(2,N) = G*WK(1,N-1) + WK(2,N)
     320  IF (WK(2,N) .EQ. ZERO)   GO TO 5008
     321  D(1,N) = (G*D(1,N-1) + D(1,N))/WK(2,N)
     322  !
     323  !  CARRY OUT BACK SUBSTITUTION
     324  !
     325   30   CONTINUE
     326  DO J=NM1,1,-1
     327     IF (WK(2,J) .EQ. ZERO)  GO TO 5008
     328     D(1,J) = (D(1,J) - WK(1,J)*D(1,J+1))/WK(2,J)
     329  END DO
     330  ! --------------------(  END  CODING FROM CUBSPL )--------------------
     331  !
     332  !  NORMAL RETURN.
     333  !
     334  RETURN
     335  !
     336  !  ERROR RETURNS.
     337  !
     338 5001   CONTINUE
     339  ! N.LT.2 RETURN.
     340  IERR = -1
     341  CALL XERMSG ('SLATEC', 'PCHSP', &
     342        'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
     343  RETURN
     344  !
     345 5002   CONTINUE
     346  ! INCFD.LT.1 RETURN.
     347  IERR = -2
     348  CALL XERMSG ('SLATEC', 'PCHSP', 'INCREMENT LESS THAN ONE', IERR, &
     349        1)
     350  RETURN
     351  !
     352 5003   CONTINUE
     353  ! X-ARRAY NOT STRICTLY INCREASING.
     354  IERR = -3
     355  CALL XERMSG ('SLATEC', 'PCHSP', 'X-ARRAY NOT STRICTLY INCREASING' &
     356        , IERR, 1)
     357  RETURN
     358  !
     359 5004   CONTINUE
     360  ! IC OUT OF RANGE RETURN.
     361  IERR = IERR - 3
     362  CALL XERMSG ('SLATEC', 'PCHSP', 'IC OUT OF RANGE', IERR, 1)
     363  RETURN
     364  !
     365 5007   CONTINUE
     366  ! NWK TOO SMALL RETURN.
     367  IERR = -7
     368  CALL XERMSG ('SLATEC', 'PCHSP', 'WORK ARRAY TOO SMALL', IERR, 1)
     369  RETURN
     370  !
     371 5008   CONTINUE
     372  ! SINGULAR SYSTEM.
     373  !   *** THEORETICALLY, THIS CAN ONLY OCCUR IF SUCCESSIVE X-VALUES   ***
     374  !   *** ARE EQUAL, WHICH SHOULD ALREADY HAVE BEEN CAUGHT (IERR=-3). ***
     375  IERR = -8
     376  CALL XERMSG ('SLATEC', 'PCHSP', 'SINGULAR LINEAR SYSTEM', IERR, &
     377        1)
     378  RETURN
     379  !
     380 5009   CONTINUE
     381  ! ERROR RETURN FROM PCHDF.
     382  !   *** THIS CASE SHOULD NEVER OCCUR ***
     383  IERR = -9
     384  CALL XERMSG ('SLATEC', 'PCHSP', 'ERROR RETURN FROM PCHDF', IERR, &
     385        1)
     386  RETURN
     387  !------------- LAST LINE OF PCHSP FOLLOWS ------------------------------
     388END SUBROUTINE PCHSP

  • LMDZ6/trunk/libf/misc/q_sat.f90

    r5245 r5246  
    55!
    66
    7       subroutine q_sat(np,temp,pres,qsat)
    8 !
    9       IMPLICIT none
    10 !======================================================================
    11 ! Autheur(s): Z.X. Li (LMD/CNRS)
    12 !  reecriture vectorisee par F. Hourdin.
    13 ! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
    14 !======================================================================
    15 ! Arguments:
    16 ! kelvin---input-R: temperature en Kelvin
    17 ! millibar--input-R: pression en mb
    18 !
    19 ! q_sat----output-R: vapeur d'eau saturante en kg/kg
    20 !======================================================================
    21 !
    22       integer np
    23       REAL temp(np),pres(np),qsat(np)
    24 !
    25       REAL r2es
    26       PARAMETER (r2es=611.14 *18.0153/28.9644)
    27 !
    28       REAL r3les, r3ies, r3es
    29       PARAMETER (R3LES=17.269)
    30       PARAMETER (R3IES=21.875)
    31 !
    32       REAL r4les, r4ies, r4es
    33       PARAMETER (R4LES=35.86)
    34       PARAMETER (R4IES=7.66)
    35 !
    36       REAL rtt
    37       PARAMETER (rtt=273.16)
    38 !
    39       REAL retv
    40       PARAMETER (retv=28.9644/18.0153 - 1.0)
     7subroutine q_sat(np,temp,pres,qsat)
     8  !
     9  IMPLICIT none
     10  !======================================================================
     11  ! Autheur(s): Z.X. Li (LMD/CNRS)
     12  !  reecriture vectorisee par F. Hourdin.
     13  ! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
     14  !======================================================================
     15  ! Arguments:
     16  ! kelvin---input-R: temperature en Kelvin
     17  ! millibar--input-R: pression en mb
     18  !
     19  ! q_sat----output-R: vapeur d'eau saturante en kg/kg
     20  !======================================================================
     21  !
     22  integer :: np
     23  REAL :: temp(np),pres(np),qsat(np)
     24  !
     25  REAL :: r2es
     26  PARAMETER (r2es=611.14 *18.0153/28.9644)
     27  !
     28  REAL :: r3les, r3ies, r3es
     29  PARAMETER (R3LES=17.269)
     30  PARAMETER (R3IES=21.875)
     31  !
     32  REAL :: r4les, r4ies, r4es
     33  PARAMETER (R4LES=35.86)
     34  PARAMETER (R4IES=7.66)
     35  !
     36  REAL :: rtt
     37  PARAMETER (rtt=273.16)
     38  !
     39  REAL :: retv
     40  PARAMETER (retv=28.9644/18.0153 - 1.0)
    4141
    42       real zqsat
    43       integer ip
    44 !
    45 !    ------------------------------------------------------------------
    46 !
    47 !
     42  real :: zqsat
     43  integer :: ip
     44  !
     45  ! ------------------------------------------------------------------
     46  !
     47  !
    4848
    49       do ip=1,np
     49  do ip=1,np
    5050
    51 !      write(*,*)'kelvin,millibar=',kelvin,millibar
    52 !       write(*,*)'temp,pres=',temp(ip),pres(ip)
    53 !
    54          IF (temp(ip) .LE. rtt) THEN
    55             r3es = r3ies
    56             r4es = r4ies
    57          ELSE
    58             r3es = r3les
    59             r4es = r4les
    60          ENDIF
    61 !
    62          zqsat=r2es/pres(ip)*EXP(r3es*(temp(ip)-rtt)/(temp(ip)-r4es))
    63          zqsat=MIN(0.5,ZQSAT)
    64          zqsat=zqsat/(1.-retv *zqsat)
    65 !
    66          qsat(ip)= zqsat
    67 !      write(*,*)'qsat=',qsat(ip)
     51   ! write(*,*)'kelvin,millibar=',kelvin,millibar
     52   !  write(*,*)'temp,pres=',temp(ip),pres(ip)
     53  !
     54     IF (temp(ip) .LE. rtt) THEN
     55        r3es = r3ies
     56        r4es = r4ies
     57     ELSE
     58        r3es = r3les
     59        r4es = r4les
     60     ENDIF
     61  !
     62     zqsat=r2es/pres(ip)*EXP(r3es*(temp(ip)-rtt)/(temp(ip)-r4es))
     63     zqsat=MIN(0.5,ZQSAT)
     64     zqsat=zqsat/(1.-retv *zqsat)
     65  !
     66     qsat(ip)= zqsat
     67   ! write(*,*)'qsat=',qsat(ip)
    6868
    69       enddo
    70 !
    71       RETURN
    72       END
     69  enddo
     70  !
     71  RETURN
     72END SUBROUTINE q_sat

  • LMDZ6/trunk/libf/misc/ran1.f90

    r5245 r5246  
    22! $Id$
    33!
    4       FUNCTION RAN1(IDUM)
    5       IMPLICIT NONE
    6       REAL RAN1
    7       REAL,SAVE :: R(97)
    8       REAL,PARAMETER :: RM1=3.8580247E-6,RM2=7.4373773E-6
    9       INTEGER,SAVE :: IFF=0
    10       integer,save :: ix1,ix2,ix3
    11       INTEGER,PARAMETER :: M1=259200,IA1=7141,IC1=54773
    12       INTEGER,PARAMETER :: M2=134456,IA2=8121,IC2=28411
    13       INTEGER,PARAMETER :: M3=243000,IA3=4561,IC3=51349
    14       INTEGER :: IDUM,J
     4FUNCTION RAN1(IDUM)
     5  IMPLICIT NONE
     6  REAL :: RAN1
     7  REAL,SAVE :: R(97)
     8  REAL,PARAMETER :: RM1=3.8580247E-6,RM2=7.4373773E-6
     9  INTEGER,SAVE :: IFF=0
     10  integer,save :: ix1,ix2,ix3
     11  INTEGER,PARAMETER :: M1=259200,IA1=7141,IC1=54773
     12  INTEGER,PARAMETER :: M2=134456,IA2=8121,IC2=28411
     13  INTEGER,PARAMETER :: M3=243000,IA3=4561,IC3=51349
     14  INTEGER :: IDUM,J
    1515
    16       IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
    17         IFF=1
    18         IX1=MOD(IC1-IDUM,M1)
    19         IX1=MOD(IA1*IX1+IC1,M1)
    20         IX2=MOD(IX1,M2)
    21         IX1=MOD(IA1*IX1+IC1,M1)
    22         IX3=MOD(IX1,M3)
    23         DO 11 J=1,97
    24           IX1=MOD(IA1*IX1+IC1,M1)
    25           IX2=MOD(IA2*IX2+IC2,M2)
    26           R(J)=(REAL(IX1)+REAL(IX2)*RM2)*RM1
    27 11      CONTINUE
    28         IDUM=1
    29       ENDIF
     16  IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
     17    IFF=1
     18    IX1=MOD(IC1-IDUM,M1)
     19    IX1=MOD(IA1*IX1+IC1,M1)
     20    IX2=MOD(IX1,M2)
     21    IX1=MOD(IA1*IX1+IC1,M1)
     22    IX3=MOD(IX1,M3)
     23    DO J=1,97
    3024      IX1=MOD(IA1*IX1+IC1,M1)
    3125      IX2=MOD(IA2*IX2+IC2,M2)
    32       IX3=MOD(IA3*IX3+IC3,M3)
    33       J=1+(97*IX3)/M3
    34       IF(J.GT.97.OR.J.LT.1) stop 1
    35       RAN1=R(J)
    3626      R(J)=(REAL(IX1)+REAL(IX2)*RM2)*RM1
    37       RETURN
    38       END
     27    END DO
     28    IDUM=1
     29  ENDIF
     30  IX1=MOD(IA1*IX1+IC1,M1)
     31  IX2=MOD(IA2*IX2+IC2,M2)
     32  IX3=MOD(IA3*IX3+IC3,M3)
     33  J=1+(97*IX3)/M3
     34  IF(J.GT.97.OR.J.LT.1) stop 1
     35  RAN1=R(J)
     36  R(J)=(REAL(IX1)+REAL(IX2)*RM2)*RM1
     37  RETURN
     38END FUNCTION RAN1

  • LMDZ6/trunk/libf/misc/sort.f90

    r5245 r5246  
    22! $Header$
    33!
    4 C
    5 C
    6       SUBROUTINE sort(n,d)
    7 c
    8 c    P.Le Van
    9 c     
    10 c...  cette routine met le tableau d  dans l'ordre croissant  ....
    11 cc   ( pour avoir l'ordre decroissant,il suffit de remplacer l'instruc
    12 c      tion  situee + bas  IF(d(j).LE.p)  THEN     par
    13 c                           IF(d(j).GE.p)  THEN
    14 c
     4!
     5!
     6SUBROUTINE sort(n,d)
     7  !
     8  ! P.Le Van
     9  !
     10  !...  cette routine met le tableau d  dans l'ordre croissant  ....
     11  !c   ( pour avoir l'ordre decroissant,il suffit de remplacer l'instruc
     12  !  tion  situee + bas  IF(d(j).LE.p)  THEN     par
     13  !                       IF(d(j).GE.p)  THEN
     14  !
    1515
    16       INTEGER n
    17       REAL d(n) , p
    18       INTEGER i,j,k
     16  INTEGER :: n
     17  REAL :: d(n) , p
     18  INTEGER :: i,j,k
    1919
    20       DO i=1,n-1
    21         k=i
    22         p=d(i)
    23         DO j=i+1,n
    24          IF(d(j).LE.p) THEN
    25            k=j
    26            p=d(j)
    27          ENDIF
    28         ENDDO
     20  DO i=1,n-1
     21    k=i
     22    p=d(i)
     23    DO j=i+1,n
     24     IF(d(j).LE.p) THEN
     25       k=j
     26       p=d(j)
     27     ENDIF
     28    ENDDO
    2929
    30        IF(k.ne.i) THEN
    31          d(k)=d(i)
    32          d(i)=p
    33        ENDIF
    34       ENDDO
     30   IF(k.ne.i) THEN
     31     d(k)=d(i)
     32     d(i)=p
     33   ENDIF
     34  ENDDO
    3535
    36        RETURN
    37        END
     36   RETURN
     37END SUBROUTINE sort

  • LMDZ6/trunk/libf/misc/xercnt.f90

    r5245 r5246  
    1 *DECK XERCNT
    2       SUBROUTINE XERCNT (LIBRAR, SUBROU, MESSG, NERR, LEVEL, KONTRL)
    3       IMPLICIT NONE
    4 C***BEGIN PROLOGUE  XERCNT
    5 C***SUBSIDIARY
    6 C***PURPOSE  Allow user control over handling of errors.
    7 C***LIBRARY   SLATEC (XERROR)
    8 C***CATEGORY  R3C
    9 C***TYPE      ALL (XERCNT-A)
    10 C***KEYWORDS  ERROR, XERROR
    11 C***AUTHOR  Jones, R. E., (SNLA)
    12 C***DESCRIPTION
    13 C
    14 C    Abstract
    15 C        Allows user control over handling of individual errors.
    16 C        Just after each message is recorded, but before it is
    17 C        processed any further (i.e., before it is printed or
    18 C        a decision to abort is made), a call is made to XERCNT.
    19 C        If the user has provided his own version of XERCNT, he
    20 C        can then override the value of KONTROL used in processing
    21 C        this message by redefining its value.
    22 C        KONTRL may be set to any value from -2 to 2.
    23 C        The meanings for KONTRL are the same as in XSETF, except
    24 C        that the value of KONTRL changes only for this message.
    25 C        If KONTRL is set to a value outside the range from -2 to 2,
    26 C        it will be moved back into that range.
    27 C
    28 C    Description of Parameters
    29 C
    30 C      --Input--
    31 C        LIBRAR - the library that the routine is in.
    32 C        SUBROU - the subroutine that XERMSG is being called from
    33 C        MESSG  - the first 20 characters of the error message.
    34 C        NERR   - same as in the call to XERMSG.
    35 C        LEVEL  - same as in the call to XERMSG.
    36 C        KONTRL - the current value of the control flag as set
    37 C                 by a call to XSETF.
    38 C
    39 C      --Output--
    40 C        KONTRL - the new value of KONTRL.  If KONTRL is not
    41 C                 defined, it will remain at its original value.
    42 C                 This changed value of control affects only
    43 C                 the current occurrence of the current message.
    44 C
    45 C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
    46 C                 Error-handling Package, SAND82-0800, Sandia
    47 C                 Laboratories, 1982.
    48 C***ROUTINES CALLED  (NONE)
    49 C***REVISION HISTORY  (YYMMDD)
    50 C   790801  DATE WRITTEN
    51 C   861211  REVISION DATE from Version 3.2
    52 C   891214  Prologue converted to Version 4.0 format.  (BAB)
    53 C   900206  Routine changed from user-callable to subsidiary.  (WRB)
    54 C   900510  Changed calling sequence to include LIBRARY and SUBROUTINE
    55 C           names, changed routine name from XERCTL to XERCNT.  (RWC)
    56 C   920501  Reformatted the REFERENCES section.  (WRB)
    57 C***END PROLOGUE  XERCNT
    58       CHARACTER*(*) LIBRAR, SUBROU, MESSG
    59       INTEGER NERR, LEVEL, KONTRL
    60 C***FIRST EXECUTABLE STATEMENT  XERCNT
    61       RETURN
    62       END
     1!DECK XERCNT
     2SUBROUTINE XERCNT (LIBRAR, SUBROU, MESSG, NERR, LEVEL, KONTRL)
     3  IMPLICIT NONE
     4  !***BEGIN PROLOGUE  XERCNT
     5  !***SUBSIDIARY
     6  !***PURPOSE  Allow user control over handling of errors.
     7  !***LIBRARY   SLATEC (XERROR)
     8  !***CATEGORY  R3C
     9  !***TYPE      ALL (XERCNT-A)
     10  !***KEYWORDS  ERROR, XERROR
     11  !***AUTHOR  Jones, R. E., (SNLA)
     12  !***DESCRIPTION
     13  !
     14  ! Abstract
     15  !    Allows user control over handling of individual errors.
     16  !    Just after each message is recorded, but before it is
     17  !    processed any further (i.e., before it is printed or
     18  !    a decision to abort is made), a call is made to XERCNT.
     19  !    If the user has provided his own version of XERCNT, he
     20  !    can then override the value of KONTROL used in processing
     21  !    this message by redefining its value.
     22  !    KONTRL may be set to any value from -2 to 2.
     23  !    The meanings for KONTRL are the same as in XSETF, except
     24  !    that the value of KONTRL changes only for this message.
     25  !    If KONTRL is set to a value outside the range from -2 to 2,
     26  !    it will be moved back into that range.
     27  !
     28  ! Description of Parameters
     29  !
     30  !  --Input--
     31  !    LIBRAR - the library that the routine is in.
     32  !    SUBROU - the subroutine that XERMSG is being called from
     33  !    MESSG  - the first 20 characters of the error message.
     34  !    NERR   - same as in the call to XERMSG.
     35  !    LEVEL  - same as in the call to XERMSG.
     36  !    KONTRL - the current value of the control flag as set
     37  !             by a call to XSETF.
     38  !
     39  !  --Output--
     40  !    KONTRL - the new value of KONTRL.  If KONTRL is not
     41  !             defined, it will remain at its original value.
     42  !             This changed value of control affects only
     43  !             the current occurrence of the current message.
     44  !
     45  !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
     46  !             Error-handling Package, SAND82-0800, Sandia
     47  !             Laboratories, 1982.
     48  !***ROUTINES CALLED  (NONE)
     49  !***REVISION HISTORY  (YYMMDD)
     50  !   790801  DATE WRITTEN
     51  !   861211  REVISION DATE from Version 3.2
     52  !   891214  Prologue converted to Version 4.0 format.  (BAB)
     53  !   900206  Routine changed from user-callable to subsidiary.  (WRB)
     54  !   900510  Changed calling sequence to include LIBRARY and SUBROUTINE
     55  !       names, changed routine name from XERCTL to XERCNT.  (RWC)
     56  !   920501  Reformatted the REFERENCES section.  (WRB)
     57  !***END PROLOGUE  XERCNT
     58  CHARACTER(len=*) :: LIBRAR, SUBROU, MESSG
     59  INTEGER :: NERR, LEVEL, KONTRL
     60  !***FIRST EXECUTABLE STATEMENT  XERCNT
     61  RETURN
     62END SUBROUTINE XERCNT

  • LMDZ6/trunk/libf/misc/xerhlt.f90

    r5245 r5246  
    1 *DECK XERHLT
    2       SUBROUTINE XERHLT (MESSG)
    3 C***BEGIN PROLOGUE  XERHLT
    4 C***SUBSIDIARY
    5 C***PURPOSE  Abort program execution and print error message.
    6 C***LIBRARY   SLATEC (XERROR)
    7 C***CATEGORY  R3C
    8 C***TYPE      ALL (XERHLT-A)
    9 C***KEYWORDS  ABORT PROGRAM EXECUTION, ERROR, XERROR
    10 C***AUTHOR  Jones, R. E., (SNLA)
    11 C***DESCRIPTION
    12 C
    13 C    Abstract
    14 C        ***Note*** machine dependent routine
    15 C        XERHLT aborts the execution of the program.
    16 C        The error message causing the abort is given in the calling
    17 C        sequence, in case one needs it for printing on a dayfile,
    18 C        for example.
    19 C
    20 C    Description of Parameters
    21 C        MESSG is as in XERMSG.
    22 C
    23 C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
    24 C                 Error-handling Package, SAND82-0800, Sandia
    25 C                 Laboratories, 1982.
    26 C***ROUTINES CALLED  (NONE)
    27 C***REVISION HISTORY  (YYMMDD)
    28 C   790801  DATE WRITTEN
    29 C   861211  REVISION DATE from Version 3.2
    30 C   891214  Prologue converted to Version 4.0 format.  (BAB)
    31 C   900206  Routine changed from user-callable to subsidiary.  (WRB)
    32 C   900510  Changed calling sequence to delete length of character
    33 C           and changed routine name from XERABT to XERHLT.  (RWC)
    34 C   920501  Reformatted the REFERENCES section.  (WRB)
    35 C***END PROLOGUE  XERHLT
    36       CHARACTER*(*) MESSG
    37 C***FIRST EXECUTABLE STATEMENT  XERHLT
    38       STOP
    39       END
     1!DECK XERHLT
     2SUBROUTINE XERHLT (MESSG)
     3  !***BEGIN PROLOGUE  XERHLT
     4  !***SUBSIDIARY
     5  !***PURPOSE  Abort program execution and print error message.
     6  !***LIBRARY   SLATEC (XERROR)
     7  !***CATEGORY  R3C
     8  !***TYPE      ALL (XERHLT-A)
     9  !***KEYWORDS  ABORT PROGRAM EXECUTION, ERROR, XERROR
     10  !***AUTHOR  Jones, R. E., (SNLA)
     11  !***DESCRIPTION
     12  !
     13  ! Abstract
     14  !    ***Note*** machine dependent routine
     15  !    XERHLT aborts the execution of the program.
     16  !    The error message causing the abort is given in the calling
     17  !    sequence, in case one needs it for printing on a dayfile,
     18  !    for example.
     19  !
     20  ! Description of Parameters
     21  !    MESSG is as in XERMSG.
     22  !
     23  !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
     24  !             Error-handling Package, SAND82-0800, Sandia
     25  !             Laboratories, 1982.
     26  !***ROUTINES CALLED  (NONE)
     27  !***REVISION HISTORY  (YYMMDD)
     28  !   790801  DATE WRITTEN
     29  !   861211  REVISION DATE from Version 3.2
     30  !   891214  Prologue converted to Version 4.0 format.  (BAB)
     31  !   900206  Routine changed from user-callable to subsidiary.  (WRB)
     32  !   900510  Changed calling sequence to delete length of character
     33  !       and changed routine name from XERABT to XERHLT.  (RWC)
     34  !   920501  Reformatted the REFERENCES section.  (WRB)
     35  !***END PROLOGUE  XERHLT
     36  CHARACTER(len=*) :: MESSG
     37  !***FIRST EXECUTABLE STATEMENT  XERHLT
     38  STOP
     39END SUBROUTINE XERHLT

  • LMDZ6/trunk/libf/misc/xermsg.f90

    r5245 r5246  
    1 *DECK XERMSG
    2       SUBROUTINE XERMSG (LIBRAR, SUBROU, MESSG, NERR, LEVEL)
    3       IMPLICIT NONE
    4 C***BEGIN PROLOGUE  XERMSG
    5 C***PURPOSE  Process error messages for SLATEC and other libraries.
    6 C***LIBRARY   SLATEC (XERROR)
    7 C***CATEGORY  R3C
    8 C***TYPE      ALL (XERMSG-A)
    9 C***KEYWORDS  ERROR MESSAGE, XERROR
    10 C***AUTHOR  Fong, Kirby, (NMFECC at LLNL)
    11 C***DESCRIPTION
    12 C
    13 C   XERMSG processes a diagnostic message in a manner determined by the
    14 C   value of LEVEL and the current value of the library error control
    15 C   flag, KONTRL.  See subroutine XSETF for details.
    16 C
    17 C    LIBRAR   A character constant (or character variable) with the name
    18 C             of the library.  This will be 'SLATEC' for the SLATEC
    19 C             Common Math Library.  The error handling package is
    20 C             general enough to be used by many libraries
    21 C             simultaneously, so it is desirable for the routine that
    22 C             detects and reports an error to identify the library name
    23 C             as well as the routine name.
    24 C
    25 C    SUBROU   A character constant (or character variable) with the name
    26 C             of the routine that detected the error.  Usually it is the
    27 C             name of the routine that is calling XERMSG.  There are
    28 C             some instances where a user callable library routine calls
    29 C             lower level subsidiary routines where the error is
    30 C             detected.  In such cases it may be more informative to
    31 C             supply the name of the routine the user called rather than
    32 C             the name of the subsidiary routine that detected the
    33 C             error.
    34 C
    35 C    MESSG    A character constant (or character variable) with the text
    36 C             of the error or warning message.  In the example below,
    37 C             the message is a character constant that contains a
    38 C             generic message.
    39 C
    40 C                   CALL XERMSG ('SLATEC', 'MMPY',
    41 C                  *'THE ORDER OF THE MATRIX EXCEEDS THE ROW DIMENSION',
    42 C                  *3, 1)
    43 C
    44 C             It is possible (and is sometimes desirable) to generate a
    45 C             specific message--e.g., one that contains actual numeric
    46 C             values.  Specific numeric values can be converted into
    47 C             character strings using formatted WRITE statements into
    48 C             character variables.  This is called standard Fortran
    49 C             internal file I/O and is exemplified in the first three
    50 C             lines of the following example.  You can also catenate
    51 C             substrings of characters to construct the error message.
    52 C             Here is an example showing the use of both writing to
    53 C             an internal file and catenating character strings.
    54 C
    55 C                   CHARACTER*5 CHARN, CHARL
    56 C                   WRITE (CHARN,10) N
    57 C                   WRITE (CHARL,10) LDA
    58 C                10 FORMAT(I5)
    59 C                   CALL XERMSG ('SLATEC', 'MMPY', 'THE ORDER'//CHARN//
    60 C                  *   ' OF THE MATRIX EXCEEDS ITS ROW DIMENSION OF'//
    61 C                  *   CHARL, 3, 1)
    62 C
    63 C             There are two subtleties worth mentioning.  One is that
    64 C             the // for character catenation is used to construct the
    65 C             error message so that no single character constant is
    66 C             continued to the next line.  This avoids confusion as to
    67 C             whether there are trailing blanks at the end of the line.
    68 C             The second is that by catenating the parts of the message
    69 C             as an actual argument rather than encoding the entire
    70 C             message into one large character variable, we avoid
    71 C             having to know how long the message will be in order to
    72 C             declare an adequate length for that large character
    73 C             variable.  XERMSG calls XERPRN to print the message using
    74 C             multiple lines if necessary.  If the message is very long,
    75 C             XERPRN will break it into pieces of 72 characters (as
    76 C             requested by XERMSG) for printing on multiple lines.
    77 C             Also, XERMSG asks XERPRN to prefix each line with ' *  '
    78 C             so that the total line length could be 76 characters.
    79 C             Note also that XERPRN scans the error message backwards
    80 C             to ignore trailing blanks.  Another feature is that
    81 C             the substring '$$' is treated as a new line sentinel
    82 C             by XERPRN.  If you want to construct a multiline
    83 C             message without having to count out multiples of 72
    84 C             characters, just use '$$' as a separator.  '$$'
    85 C             obviously must occur within 72 characters of the
    86 C             start of each line to have its intended effect since
    87 C             XERPRN is asked to wrap around at 72 characters in
    88 C             addition to looking for '$$'.
    89 C
    90 C    NERR     An integer value that is chosen by the library routine's
    91 C             author.  It must be in the range -99 to 999 (three
    92 C             printable digits).  Each distinct error should have its
    93 C             own error number.  These error numbers should be described
    94 C             in the machine readable documentation for the routine.
    95 C             The error numbers need be unique only within each routine,
    96 C             so it is reasonable for each routine to start enumerating
    97 C             errors from 1 and proceeding to the next integer.
    98 C
    99 C    LEVEL    An integer value in the range 0 to 2 that indicates the
    100 C             level (severity) of the error.  Their meanings are
    101 C
    102 C            -1  A warning message.  This is used if it is not clear
    103 C                that there really is an error, but the user's attention
    104 C                may be needed.  An attempt is made to only print this
    105 C                message once.
    106 C
    107 C             0  A warning message.  This is used if it is not clear
    108 C                that there really is an error, but the user's attention
    109 C                may be needed.
    110 C
    111 C             1  A recoverable error.  This is used even if the error is
    112 C                so serious that the routine cannot return any useful
    113 C                answer.  If the user has told the error package to
    114 C                return after recoverable errors, then XERMSG will
    115 C                return to the Library routine which can then return to
    116 C                the user's routine.  The user may also permit the error
    117 C                package to terminate the program upon encountering a
    118 C                recoverable error.
    119 C
    120 C             2  A fatal error.  XERMSG will not return to its caller
    121 C                after it receives a fatal error.  This level should
    122 C                hardly ever be used; it is much better to allow the
    123 C                user a chance to recover.  An example of one of the few
    124 C                cases in which it is permissible to declare a level 2
    125 C                error is a reverse communication Library routine that
    126 C                is likely to be called repeatedly until it integrates
    127 C                across some interval.  If there is a serious error in
    128 C                the input such that another step cannot be taken and
    129 C                the Library routine is called again without the input
    130 C                error having been corrected by the caller, the Library
    131 C                routine will probably be called forever with improper
    132 C                input.  In this case, it is reasonable to declare the
    133 C                error to be fatal.
    134 C
    135 C    Each of the arguments to XERMSG is input; none will be modified by
    136 C    XERMSG.  A routine may make multiple calls to XERMSG with warning
    137 C    level messages; however, after a call to XERMSG with a recoverable
    138 C    error, the routine should return to the user.  Do not try to call
    139 C    XERMSG with a second recoverable error after the first recoverable
    140 C    error because the error package saves the error number.  The user
    141 C    can retrieve this error number by calling another entry point in
    142 C    the error handling package and then clear the error number when
    143 C    recovering from the error.  Calling XERMSG in succession causes the
    144 C    old error number to be overwritten by the latest error number.
    145 C    This is considered harmless for error numbers associated with
    146 C    warning messages but must not be done for error numbers of serious
    147 C    errors.  After a call to XERMSG with a recoverable error, the user
    148 C    must be given a chance to call NUMXER or XERCLR to retrieve or
    149 C    clear the error number.
    150 C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
    151 C                 Error-handling Package, SAND82-0800, Sandia
    152 C                 Laboratories, 1982.
    153 C***ROUTINES CALLED  FDUMP, J4SAVE, XERCNT, XERHLT, XERPRN, XERSVE
    154 C***REVISION HISTORY  (YYMMDD)
    155 C   880101  DATE WRITTEN
    156 C   880621  REVISED AS DIRECTED AT SLATEC CML MEETING OF FEBRUARY 1988.
    157 C           THERE ARE TWO BASIC CHANGES.
    158 C           1.  A NEW ROUTINE, XERPRN, IS USED INSTEAD OF XERPRT TO
    159 C               PRINT MESSAGES.  THIS ROUTINE WILL BREAK LONG MESSAGES
    160 C               INTO PIECES FOR PRINTING ON MULTIPLE LINES.  '$$' IS
    161 C               ACCEPTED AS A NEW LINE SENTINEL.  A PREFIX CAN BE
    162 C               ADDED TO EACH LINE TO BE PRINTED.  XERMSG USES EITHER
    163 C               ' ***' OR ' *  ' AND LONG MESSAGES ARE BROKEN EVERY
    164 C               72 CHARACTERS (AT MOST) SO THAT THE MAXIMUM LINE
    165 C               LENGTH OUTPUT CAN NOW BE AS GREAT AS 76.
    166 C           2.  THE TEXT OF ALL MESSAGES IS NOW IN UPPER CASE SINCE THE
    167 C               FORTRAN STANDARD DOCUMENT DOES NOT ADMIT THE EXISTENCE
    168 C               OF LOWER CASE.
    169 C   880708  REVISED AFTER THE SLATEC CML MEETING OF JUNE 29 AND 30.
    170 C           THE PRINCIPAL CHANGES ARE
    171 C           1.  CLARIFY COMMENTS IN THE PROLOGUES
    172 C           2.  RENAME XRPRNT TO XERPRN
    173 C           3.  REWORK HANDLING OF '$$' IN XERPRN TO HANDLE BLANK LINES
    174 C               SIMILAR TO THE WAY FORMAT STATEMENTS HANDLE THE /
    175 C               CHARACTER FOR NEW RECORDS.
    176 C   890706  REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO
    177 C           CLEAN UP THE CODING.
    178 C   890721  REVISED TO USE NEW FEATURE IN XERPRN TO COUNT CHARACTERS IN
    179 C           PREFIX.
    180 C   891013  REVISED TO CORRECT COMMENTS.
    181 C   891214  Prologue converted to Version 4.0 format.  (WRB)
    182 C   900510  Changed test on NERR to be -9999999 < NERR < 99999999, but
    183 C           NERR .ne. 0, and on LEVEL to be -2 < LEVEL < 3.  Added
    184 C           LEVEL=-1 logic, changed calls to XERSAV to XERSVE, and
    185 C           XERCTL to XERCNT.  (RWC)
    186 C   920501  Reformatted the REFERENCES section.  (WRB)
    187 C***END PROLOGUE  XERMSG
    188       CHARACTER*(*) LIBRAR, SUBROU, MESSG
    189       CHARACTER*8 XLIBR, XSUBR
    190       CHARACTER*72 TEMP
    191       CHARACTER*20 LFIRST
    192       INTEGER NERR, LEVEL, LKNTRL
    193       INTEGER J4SAVE, MAXMES, KDUMMY, I, KOUNT, LERR, LLEVEL
    194       INTEGER MKNTRL, LTEMP
    195 C***FIRST EXECUTABLE STATEMENT  XERMSG
    196       LKNTRL = J4SAVE (2, 0, .FALSE.)
    197       MAXMES = J4SAVE (4, 0, .FALSE.)
    198 C
    199 C       LKNTRL IS A LOCAL COPY OF THE CONTROL FLAG KONTRL.
    200 C       MAXMES IS THE MAXIMUM NUMBER OF TIMES ANY PARTICULAR MESSAGE
    201 C          SHOULD BE PRINTED.
    202 C
    203 C       WE PRINT A FATAL ERROR MESSAGE AND TERMINATE FOR AN ERROR IN
    204 C          CALLING XERMSG.  THE ERROR NUMBER SHOULD BE POSITIVE,
    205 C          AND THE LEVEL SHOULD BE BETWEEN 0 AND 2.
    206 C
    207       IF (NERR.LT.-9999999 .OR. NERR.GT.99999999 .OR. NERR.EQ.0 .OR.
    208      *   LEVEL.LT.-1 .OR. LEVEL.GT.2) THEN
    209          CALL XERPRN (' ***', -1, 'FATAL ERROR IN...$$ ' //
    210      *      'XERMSG -- INVALID ERROR NUMBER OR LEVEL$$ '//
    211      *      'JOB ABORT DUE TO FATAL ERROR.', 72)
    212          CALL XERSVE (' ', ' ', ' ', 0, 0, 0, KDUMMY)
    213          CALL XERHLT (' ***XERMSG -- INVALID INPUT')
    214          RETURN
    215       ENDIF
    216 C
    217 C       RECORD THE MESSAGE.
    218 C
    219       I = J4SAVE (1, NERR, .TRUE.)
    220       CALL XERSVE (LIBRAR, SUBROU, MESSG, 1, NERR, LEVEL, KOUNT)
    221 C
    222 C       HANDLE PRINT-ONCE WARNING MESSAGES.
    223 C
    224       IF (LEVEL.EQ.-1 .AND. KOUNT.GT.1) RETURN
    225 C
    226 C       ALLOW TEMPORARY USER OVERRIDE OF THE CONTROL FLAG.
    227 C
    228       XLIBR  = LIBRAR
    229       XSUBR  = SUBROU
    230       LFIRST = MESSG
    231       LERR   = NERR
    232       LLEVEL = LEVEL
    233       CALL XERCNT (XLIBR, XSUBR, LFIRST, LERR, LLEVEL, LKNTRL)
    234 C
    235       LKNTRL = MAX(-2, MIN(2,LKNTRL))
    236       MKNTRL = ABS(LKNTRL)
    237 C
    238 C       SKIP PRINTING IF THE CONTROL FLAG VALUE AS RESET IN XERCNT IS
    239 C       ZERO AND THE ERROR IS NOT FATAL.
    240 C
    241       IF (LEVEL.LT.2 .AND. LKNTRL.EQ.0) GO TO 30
    242       IF (LEVEL.EQ.0 .AND. KOUNT.GT.MAXMES) GO TO 30
    243       IF (LEVEL.EQ.1 .AND. KOUNT.GT.MAXMES .AND. MKNTRL.EQ.1) GO TO 30
    244       IF (LEVEL.EQ.2 .AND. KOUNT.GT.MAX(1,MAXMES)) GO TO 30
    245 C
    246 C       ANNOUNCE THE NAMES OF THE LIBRARY AND SUBROUTINE BY BUILDING A
    247 C       MESSAGE IN CHARACTER VARIABLE TEMP (NOT EXCEEDING 66 CHARACTERS)
    248 C       AND SENDING IT OUT VIA XERPRN.  PRINT ONLY IF CONTROL FLAG
    249 C       IS NOT ZERO.
    250 C
    251       IF (LKNTRL .NE. 0) THEN
    252          TEMP(1:21) = 'MESSAGE FROM ROUTINE '
    253          I = MIN(LEN(SUBROU), 16)
    254          TEMP(22:21+I) = SUBROU(1:I)
    255          TEMP(22+I:33+I) = ' IN LIBRARY '
    256          LTEMP = 33 + I
    257          I = MIN(LEN(LIBRAR), 16)
    258          TEMP(LTEMP+1:LTEMP+I) = LIBRAR (1:I)
    259          TEMP(LTEMP+I+1:LTEMP+I+1) = '.'
    260          LTEMP = LTEMP + I + 1
    261          CALL XERPRN (' ***', -1, TEMP(1:LTEMP), 72)
    262       ENDIF
    263 C
    264 C       IF LKNTRL IS POSITIVE, PRINT AN INTRODUCTORY LINE BEFORE
    265 C       PRINTING THE MESSAGE.  THE INTRODUCTORY LINE TELLS THE CHOICE
    266 C       FROM EACH OF THE FOLLOWING THREE OPTIONS.
    267 C       1.  LEVEL OF THE MESSAGE
    268 C              'INFORMATIVE MESSAGE'
    269 C              'POTENTIALLY RECOVERABLE ERROR'
    270 C              'FATAL ERROR'
    271 C       2.  WHETHER CONTROL FLAG WILL ALLOW PROGRAM TO CONTINUE
    272 C              'PROG CONTINUES'
    273 C              'PROG ABORTED'
    274 C       3.  WHETHER OR NOT A TRACEBACK WAS REQUESTED.  (THE TRACEBACK
    275 C           MAY NOT BE IMPLEMENTED AT SOME SITES, SO THIS ONLY TELLS
    276 C           WHAT WAS REQUESTED, NOT WHAT WAS DELIVERED.)
    277 C              'TRACEBACK REQUESTED'
    278 C              'TRACEBACK NOT REQUESTED'
    279 C       NOTICE THAT THE LINE INCLUDING FOUR PREFIX CHARACTERS WILL NOT
    280 C       EXCEED 74 CHARACTERS.
    281 C       WE SKIP THE NEXT BLOCK IF THE INTRODUCTORY LINE IS NOT NEEDED.
    282 C
    283       IF (LKNTRL .GT. 0) THEN
    284 C
    285 C       THE FIRST PART OF THE MESSAGE TELLS ABOUT THE LEVEL.
    286 C
    287          IF (LEVEL .LE. 0) THEN
    288             TEMP(1:20) = 'INFORMATIVE MESSAGE,'
    289             LTEMP = 20
    290          ELSEIF (LEVEL .EQ. 1) THEN
    291             TEMP(1:30) = 'POTENTIALLY RECOVERABLE ERROR,'
    292             LTEMP = 30
    293          ELSE
    294             TEMP(1:12) = 'FATAL ERROR,'
    295             LTEMP = 12
    296          ENDIF
    297 C
    298 C       THEN WHETHER THE PROGRAM WILL CONTINUE.
    299 C
    300          IF ((MKNTRL.EQ.2 .AND. LEVEL.GE.1) .OR.
    301      *       (MKNTRL.EQ.1 .AND. LEVEL.EQ.2)) THEN
    302             TEMP(LTEMP+1:LTEMP+14) = ' PROG ABORTED,'
    303             LTEMP = LTEMP + 14
    304          ELSE
    305             TEMP(LTEMP+1:LTEMP+16) = ' PROG CONTINUES,'
    306             LTEMP = LTEMP + 16
    307          ENDIF
    308 C
    309 C       FINALLY TELL WHETHER THERE SHOULD BE A TRACEBACK.
    310 C
    311          IF (LKNTRL .GT. 0) THEN
    312             TEMP(LTEMP+1:LTEMP+20) = ' TRACEBACK REQUESTED'
    313             LTEMP = LTEMP + 20
    314          ELSE
    315             TEMP(LTEMP+1:LTEMP+24) = ' TRACEBACK NOT REQUESTED'
    316             LTEMP = LTEMP + 24
    317          ENDIF
    318          CALL XERPRN (' ***', -1, TEMP(1:LTEMP), 72)
    319       ENDIF
    320 C
    321 C       NOW SEND OUT THE MESSAGE.
    322 C
    323       CALL XERPRN (' *  ', -1, MESSG, 72)
    324 C
    325 C       IF LKNTRL IS POSITIVE, WRITE THE ERROR NUMBER AND REQUEST A
    326 C          TRACEBACK.
    327 C
    328       IF (LKNTRL .GT. 0) THEN
    329          WRITE (TEMP, '(''ERROR NUMBER = '', I8)') NERR
    330          DO 10 I=16,22
    331             IF (TEMP(I:I) .NE. ' ') GO TO 20
    332    10    CONTINUE
    333 C
    334    20    CALL XERPRN (' *  ', -1, TEMP(1:15) // TEMP(I:23), 72)
    335          CALL FDUMP
    336       ENDIF
    337 C
    338 C       IF LKNTRL IS NOT ZERO, PRINT A BLANK LINE AND AN END OF MESSAGE.
    339 C
    340       IF (LKNTRL .NE. 0) THEN
    341          CALL XERPRN (' *  ', -1, ' ', 72)
    342          CALL XERPRN (' ***', -1, 'END OF MESSAGE', 72)
    343          CALL XERPRN ('    ',  0, ' ', 72)
    344       ENDIF
    345 C
    346 C       IF THE ERROR IS NOT FATAL OR THE ERROR IS RECOVERABLE AND THE
    347 C       CONTROL FLAG IS SET FOR RECOVERY, THEN RETURN.
    348 C
    349    30 IF (LEVEL.LE.0 .OR. (LEVEL.EQ.1 .AND. MKNTRL.LE.1)) RETURN
    350 C
    351 C       THE PROGRAM WILL BE STOPPED DUE TO AN UNRECOVERED ERROR OR A
    352 C       FATAL ERROR.  PRINT THE REASON FOR THE ABORT AND THE ERROR
    353 C       SUMMARY IF THE CONTROL FLAG AND THE MAXIMUM ERROR COUNT PERMIT.
    354 C
    355       IF (LKNTRL.GT.0 .AND. KOUNT.LT.MAX(1,MAXMES)) THEN
    356          IF (LEVEL .EQ. 1) THEN
    357             CALL XERPRN
    358      *         (' ***', -1, 'JOB ABORT DUE TO UNRECOVERED ERROR.', 72)
    359          ELSE
    360             CALL XERPRN(' ***', -1, 'JOB ABORT DUE TO FATAL ERROR.', 72)
    361          ENDIF
    362          CALL XERSVE (' ', ' ', ' ', -1, 0, 0, KDUMMY)
    363          CALL XERHLT (' ')
    364       ELSE
    365          CALL XERHLT (MESSG)
    366       ENDIF
    367       RETURN
    368       END
     1!DECK XERMSG
     2SUBROUTINE XERMSG (LIBRAR, SUBROU, MESSG, NERR, LEVEL)
     3  IMPLICIT NONE
     4  !***BEGIN PROLOGUE  XERMSG
     5  !***PURPOSE  Process error messages for SLATEC and other libraries.
     6  !***LIBRARY   SLATEC (XERROR)
     7  !***CATEGORY  R3C
     8  !***TYPE      ALL (XERMSG-A)
     9  !***KEYWORDS  ERROR MESSAGE, XERROR
     10  !***AUTHOR  Fong, Kirby, (NMFECC at LLNL)
     11  !***DESCRIPTION
     12  !
     13  !   XERMSG processes a diagnostic message in a manner determined by the
     14  !   value of LEVEL and the current value of the library error control
     15  !   flag, KONTRL.  See subroutine XSETF for details.
     16  !
     17  !    LIBRAR   A character constant (or character variable) with the name
     18  !         of the library.  This will be 'SLATEC' for the SLATEC
     19  !         Common Math Library.  The error handling package is
     20  !         general enough to be used by many libraries
     21  !         simultaneously, so it is desirable for the routine that
     22  !         detects and reports an error to identify the library name
     23  !         as well as the routine name.
     24  !
     25  !    SUBROU   A character constant (or character variable) with the name
     26  !             of the routine that detected the error.  Usually it is the
     27  !         name of the routine that is calling XERMSG.  There are
     28  !         some instances where a user callable library routine calls
     29  !         lower level subsidiary routines where the error is
     30  !         detected.  In such cases it may be more informative to
     31  !         supply the name of the routine the user called rather than
     32  !         the name of the subsidiary routine that detected the
     33  !         error.
     34  !
     35  !    MESSG    A character constant (or character variable) with the text
     36  !         of the error or warning message.  In the example below,
     37  !         the message is a character constant that contains a
     38  !         generic message.
     39  !
     40  !               CALL XERMSG ('SLATEC', 'MMPY',
     41  !              *'THE ORDER OF THE MATRIX EXCEEDS THE ROW DIMENSION',
     42  !              *3, 1)
     43  !
     44  !         It is possible (and is sometimes desirable) to generate a
     45  !         specific message--e.g., one that contains actual numeric
     46  !         values.  Specific numeric values can be converted into
     47  !         character strings using formatted WRITE statements into
     48  !         character variables.  This is called standard Fortran
     49  !         internal file I/O and is exemplified in the first three
     50  !         lines of the following example.  You can also catenate
     51  !         substrings of characters to construct the error message.
     52  !         Here is an example showing the use of both writing to
     53  !         an internal file and catenating character strings.
     54  !
     55  !               CHARACTER*5 CHARN, CHARL
     56  !               WRITE (CHARN,10) N
     57  !               WRITE (CHARL,10) LDA
     58  !            10 FORMAT(I5)
     59  !               CALL XERMSG ('SLATEC', 'MMPY', 'THE ORDER'//CHARN//
     60  !              *   ' OF THE MATRIX EXCEEDS ITS ROW DIMENSION OF'//
     61  !              *   CHARL, 3, 1)
     62  !
     63  !         There are two subtleties worth mentioning.  One is that
     64  !         the // for character catenation is used to construct the
     65  !         error message so that no single character constant is
     66  !         continued to the next line.  This avoids confusion as to
     67  !         whether there are trailing blanks at the end of the line.
     68  !         The second is that by catenating the parts of the message
     69  !             as an actual argument rather than encoding the entire
     70  !         message into one large character variable, we avoid
     71  !         having to know how long the message will be in order to
     72  !         declare an adequate length for that large character
     73  !         variable.  XERMSG calls XERPRN to print the message using
     74  !         multiple lines if necessary.  If the message is very long,
     75  !         XERPRN will break it into pieces of 72 characters (as
     76  !         requested by XERMSG) for printing on multiple lines.
     77  !         Also, XERMSG asks XERPRN to prefix each line with ' *  '
     78  !         so that the total line length could be 76 characters.
     79  !         Note also that XERPRN scans the error message backwards
     80  !         to ignore trailing blanks.  Another feature is that
     81  !         the substring '$$' is treated as a new line sentinel
     82  !         by XERPRN.  If you want to construct a multiline
     83  !         message without having to count out multiples of 72
     84  !         characters, just use '$$' as a separator.  '$$'
     85  !         obviously must occur within 72 characters of the
     86  !         start of each line to have its intended effect since
     87  !         XERPRN is asked to wrap around at 72 characters in
     88  !         addition to looking for '$$'.
     89  !
     90  !    NERR     An integer value that is chosen by the library routine's
     91  !         author.  It must be in the range -99 to 999 (three
     92  !         printable digits).  Each distinct error should have its
     93  !         own error number.  These error numbers should be described
     94  !         in the machine readable documentation for the routine.
     95  !         The error numbers need be unique only within each routine,
     96  !         so it is reasonable for each routine to start enumerating
     97  !         errors from 1 and proceeding to the next integer.
     98  !
     99  !    LEVEL    An integer value in the range 0 to 2 that indicates the
     100  !         level (severity) of the error.  Their meanings are
     101  !
     102  !        -1  A warning message.  This is used if it is not clear
     103  !            that there really is an error, but the user's attention
     104  !            may be needed.  An attempt is made to only print this
     105  !            message once.
     106  !
     107  !         0  A warning message.  This is used if it is not clear
     108  !            that there really is an error, but the user's attention
     109  !            may be needed.
     110  !
     111  !         1  A recoverable error.  This is used even if the error is
     112  !            so serious that the routine cannot return any useful
     113  !            answer.  If the user has told the error package to
     114  !            return after recoverable errors, then XERMSG will
     115  !            return to the Library routine which can then return to
     116  !            the user's routine.  The user may also permit the error
     117  !            package to terminate the program upon encountering a
     118  !            recoverable error.
     119  !
     120  !         2  A fatal error.  XERMSG will not return to its caller
     121  !            after it receives a fatal error.  This level should
     122  !            hardly ever be used; it is much better to allow the
     123  !            user a chance to recover.  An example of one of the few
     124  !            cases in which it is permissible to declare a level 2
     125  !            error is a reverse communication Library routine that
     126  !            is likely to be called repeatedly until it integrates
     127  !            across some interval.  If there is a serious error in
     128  !            the input such that another step cannot be taken and
     129  !            the Library routine is called again without the input
     130  !            error having been corrected by the caller, the Library
     131  !            routine will probably be called forever with improper
     132  !            input.  In this case, it is reasonable to declare the
     133  !            error to be fatal.
     134  !
     135  !    Each of the arguments to XERMSG is input; none will be modified by
     136  !    XERMSG.  A routine may make multiple calls to XERMSG with warning
     137  !    level messages; however, after a call to XERMSG with a recoverable
     138  !    error, the routine should return to the user.  Do not try to call
     139  !    XERMSG with a second recoverable error after the first recoverable
     140  !    error because the error package saves the error number.  The user
     141  !    can retrieve this error number by calling another entry point in
     142  !    the error handling package and then clear the error number when
     143  !    recovering from the error.  Calling XERMSG in succession causes the
     144  !    old error number to be overwritten by the latest error number.
     145  !    This is considered harmless for error numbers associated with
     146  !    warning messages but must not be done for error numbers of serious
     147  !    errors.  After a call to XERMSG with a recoverable error, the user
     148  !    must be given a chance to call NUMXER or XERCLR to retrieve or
     149  !    clear the error number.
     150  !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
     151  !             Error-handling Package, SAND82-0800, Sandia
     152  !             Laboratories, 1982.
     153  !***ROUTINES CALLED  FDUMP, J4SAVE, XERCNT, XERHLT, XERPRN, XERSVE
     154  !***REVISION HISTORY  (YYMMDD)
     155  !   880101  DATE WRITTEN
     156  !   880621  REVISED AS DIRECTED AT SLATEC CML MEETING OF FEBRUARY 1988.
     157  !       THERE ARE TWO BASIC CHANGES.
     158  !       1.  A NEW ROUTINE, XERPRN, IS USED INSTEAD OF XERPRT TO
     159  !           PRINT MESSAGES.  THIS ROUTINE WILL BREAK LONG MESSAGES
     160  !           INTO PIECES FOR PRINTING ON MULTIPLE LINES.  '$$' IS
     161  !           ACCEPTED AS A NEW LINE SENTINEL.  A PREFIX CAN BE
     162  !           ADDED TO EACH LINE TO BE PRINTED.  XERMSG USES EITHER
     163  !           ' ***' OR ' *  ' AND LONG MESSAGES ARE BROKEN EVERY
     164  !           72 CHARACTERS (AT MOST) SO THAT THE MAXIMUM LINE
     165  !           LENGTH OUTPUT CAN NOW BE AS GREAT AS 76.
     166  !       2.  THE TEXT OF ALL MESSAGES IS NOW IN UPPER CASE SINCE THE
     167  !           FORTRAN STANDARD DOCUMENT DOES NOT ADMIT THE EXISTENCE
     168  !           OF LOWER CASE.
     169  !   880708  REVISED AFTER THE SLATEC CML MEETING OF JUNE 29 AND 30.
     170  !       THE PRINCIPAL CHANGES ARE
     171  !       1.  CLARIFY COMMENTS IN THE PROLOGUES
     172  !       2.  RENAME XRPRNT TO XERPRN
     173  !       3.  REWORK HANDLING OF '$$' IN XERPRN TO HANDLE BLANK LINES
     174  !           SIMILAR TO THE WAY FORMAT STATEMENTS HANDLE THE /
     175  !           CHARACTER FOR NEW RECORDS.
     176  !   890706  REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO
     177  !       CLEAN UP THE CODING.
     178  !   890721  REVISED TO USE NEW FEATURE IN XERPRN TO COUNT CHARACTERS IN
     179  !       PREFIX.
     180  !   891013  REVISED TO CORRECT COMMENTS.
     181  !   891214  Prologue converted to Version 4.0 format.  (WRB)
     182  !   900510  Changed test on NERR to be -9999999 < NERR < 99999999, but
     183  !       NERR .ne. 0, and on LEVEL to be -2 < LEVEL < 3.  Added
     184  !       LEVEL=-1 logic, changed calls to XERSAV to XERSVE, and
     185  !       XERCTL to XERCNT.  (RWC)
     186  !   920501  Reformatted the REFERENCES section.  (WRB)
     187  !***END PROLOGUE  XERMSG
     188  CHARACTER(len=*) :: LIBRAR, SUBROU, MESSG
     189  CHARACTER(len=8) :: XLIBR, XSUBR
     190  CHARACTER(len=72) :: TEMP
     191  CHARACTER(len=20) :: LFIRST
     192  INTEGER :: NERR, LEVEL, LKNTRL
     193  INTEGER :: J4SAVE, MAXMES, KDUMMY, I, KOUNT, LERR, LLEVEL
     194  INTEGER :: MKNTRL, LTEMP
     195  !***FIRST EXECUTABLE STATEMENT  XERMSG
     196  LKNTRL = J4SAVE (2, 0, .FALSE.)
     197  MAXMES = J4SAVE (4, 0, .FALSE.)
     198  !
     199  !   LKNTRL IS A LOCAL COPY OF THE CONTROL FLAG KONTRL.
     200  !   MAXMES IS THE MAXIMUM NUMBER OF TIMES ANY PARTICULAR MESSAGE
     201  !      SHOULD BE PRINTED.
     202  !
     203  !   WE PRINT A FATAL ERROR MESSAGE AND TERMINATE FOR AN ERROR IN
     204  !      CALLING XERMSG.  THE ERROR NUMBER SHOULD BE POSITIVE,
     205  !      AND THE LEVEL SHOULD BE BETWEEN 0 AND 2.
     206  !
     207  IF (NERR.LT.-9999999 .OR. NERR.GT.99999999 .OR. NERR.EQ.0 .OR. &
     208        LEVEL.LT.-1 .OR. LEVEL.GT.2) THEN
     209     CALL XERPRN (' ***', -1, 'FATAL ERROR IN...$$ ' // &
     210           'XERMSG -- INVALID ERROR NUMBER OR LEVEL$$ '// &
     211           'JOB ABORT DUE TO FATAL ERROR.', 72)
     212     CALL XERSVE (' ', ' ', ' ', 0, 0, 0, KDUMMY)
     213     CALL XERHLT (' ***XERMSG -- INVALID INPUT')
     214     RETURN
     215  ENDIF
     216  !
     217  !   RECORD THE MESSAGE.
     218  !
     219  I = J4SAVE (1, NERR, .TRUE.)
     220  CALL XERSVE (LIBRAR, SUBROU, MESSG, 1, NERR, LEVEL, KOUNT)
     221  !
     222  !   HANDLE PRINT-ONCE WARNING MESSAGES.
     223  !
     224  IF (LEVEL.EQ.-1 .AND. KOUNT.GT.1) RETURN
     225  !
     226  !   ALLOW TEMPORARY USER OVERRIDE OF THE CONTROL FLAG.
     227  !
     228  XLIBR  = LIBRAR
     229  XSUBR  = SUBROU
     230  LFIRST = MESSG
     231  LERR   = NERR
     232  LLEVEL = LEVEL
     233  CALL XERCNT (XLIBR, XSUBR, LFIRST, LERR, LLEVEL, LKNTRL)
     234  !
     235  LKNTRL = MAX(-2, MIN(2,LKNTRL))
     236  MKNTRL = ABS(LKNTRL)
     237  !
     238  !   SKIP PRINTING IF THE CONTROL FLAG VALUE AS RESET IN XERCNT IS
     239  !   ZERO AND THE ERROR IS NOT FATAL.
     240  !
     241  IF (LEVEL.LT.2 .AND. LKNTRL.EQ.0) GO TO 30
     242  IF (LEVEL.EQ.0 .AND. KOUNT.GT.MAXMES) GO TO 30
     243  IF (LEVEL.EQ.1 .AND. KOUNT.GT.MAXMES .AND. MKNTRL.EQ.1) GO TO 30
     244  IF (LEVEL.EQ.2 .AND. KOUNT.GT.MAX(1,MAXMES)) GO TO 30
     245  !
     246  !   ANNOUNCE THE NAMES OF THE LIBRARY AND SUBROUTINE BY BUILDING A
     247  !   MESSAGE IN CHARACTER VARIABLE TEMP (NOT EXCEEDING 66 CHARACTERS)
     248  !   AND SENDING IT OUT VIA XERPRN.  PRINT ONLY IF CONTROL FLAG
     249  !   IS NOT ZERO.
     250  !
     251  IF (LKNTRL .NE. 0) THEN
     252     TEMP(1:21) = 'MESSAGE FROM ROUTINE '
     253     I = MIN(LEN(SUBROU), 16)
     254     TEMP(22:21+I) = SUBROU(1:I)
     255     TEMP(22+I:33+I) = ' IN LIBRARY '
     256     LTEMP = 33 + I
     257     I = MIN(LEN(LIBRAR), 16)
     258     TEMP(LTEMP+1:LTEMP+I) = LIBRAR (1:I)
     259     TEMP(LTEMP+I+1:LTEMP+I+1) = '.'
     260     LTEMP = LTEMP + I + 1
     261     CALL XERPRN (' ***', -1, TEMP(1:LTEMP), 72)
     262  ENDIF
     263  !
     264  !   IF LKNTRL IS POSITIVE, PRINT AN INTRODUCTORY LINE BEFORE
     265  !   PRINTING THE MESSAGE.  THE INTRODUCTORY LINE TELLS THE CHOICE
     266  !   FROM EACH OF THE FOLLOWING THREE OPTIONS.
     267  !   1.  LEVEL OF THE MESSAGE
     268  !          'INFORMATIVE MESSAGE'
     269  !          'POTENTIALLY RECOVERABLE ERROR'
     270  !          'FATAL ERROR'
     271  !   2.  WHETHER CONTROL FLAG WILL ALLOW PROGRAM TO CONTINUE
     272  !          'PROG CONTINUES'
     273  !          'PROG ABORTED'
     274  !   3.  WHETHER OR NOT A TRACEBACK WAS REQUESTED.  (THE TRACEBACK
     275  !       MAY NOT BE IMPLEMENTED AT SOME SITES, SO THIS ONLY TELLS
     276  !       WHAT WAS REQUESTED, NOT WHAT WAS DELIVERED.)
     277  !          'TRACEBACK REQUESTED'
     278  !          'TRACEBACK NOT REQUESTED'
     279  !   NOTICE THAT THE LINE INCLUDING FOUR PREFIX CHARACTERS WILL NOT
     280  !   EXCEED 74 CHARACTERS.
     281  !   WE SKIP THE NEXT BLOCK IF THE INTRODUCTORY LINE IS NOT NEEDED.
     282  !
     283  IF (LKNTRL .GT. 0) THEN
     284  !
     285  !   THE FIRST PART OF THE MESSAGE TELLS ABOUT THE LEVEL.
     286  !
     287     IF (LEVEL .LE. 0) THEN
     288        TEMP(1:20) = 'INFORMATIVE MESSAGE,'
     289        LTEMP = 20
     290     ELSEIF (LEVEL .EQ. 1) THEN
     291        TEMP(1:30) = 'POTENTIALLY RECOVERABLE ERROR,'
     292        LTEMP = 30
     293     ELSE
     294        TEMP(1:12) = 'FATAL ERROR,'
     295        LTEMP = 12
     296     ENDIF
     297  !
     298  !   THEN WHETHER THE PROGRAM WILL CONTINUE.
     299  !
     300     IF ((MKNTRL.EQ.2 .AND. LEVEL.GE.1) .OR. &
     301           (MKNTRL.EQ.1 .AND. LEVEL.EQ.2)) THEN
     302        TEMP(LTEMP+1:LTEMP+14) = ' PROG ABORTED,'
     303        LTEMP = LTEMP + 14
     304     ELSE
     305        TEMP(LTEMP+1:LTEMP+16) = ' PROG CONTINUES,'
     306        LTEMP = LTEMP + 16
     307     ENDIF
     308  !
     309  !   FINALLY TELL WHETHER THERE SHOULD BE A TRACEBACK.
     310  !
     311     IF (LKNTRL .GT. 0) THEN
     312        TEMP(LTEMP+1:LTEMP+20) = ' TRACEBACK REQUESTED'
     313        LTEMP = LTEMP + 20
     314     ELSE
     315        TEMP(LTEMP+1:LTEMP+24) = ' TRACEBACK NOT REQUESTED'
     316        LTEMP = LTEMP + 24
     317     ENDIF
     318     CALL XERPRN (' ***', -1, TEMP(1:LTEMP), 72)
     319  ENDIF
     320  !
     321  !   NOW SEND OUT THE MESSAGE.
     322  !
     323  CALL XERPRN (' *  ', -1, MESSG, 72)
     324  !
     325  !   IF LKNTRL IS POSITIVE, WRITE THE ERROR NUMBER AND REQUEST A
     326  !      TRACEBACK.
     327  !
     328  IF (LKNTRL .GT. 0) THEN
     329     WRITE (TEMP, '(''ERROR NUMBER = '', I8)') NERR
     330     DO I=16,22
     331        IF (TEMP(I:I) .NE. ' ') GO TO 20
     332     END DO
     333  !
     334   20   CALL XERPRN (' *  ', -1, TEMP(1:15) // TEMP(I:23), 72)
     335     CALL FDUMP
     336  ENDIF
     337  !
     338  !   IF LKNTRL IS NOT ZERO, PRINT A BLANK LINE AND AN END OF MESSAGE.
     339  !
     340  IF (LKNTRL .NE. 0) THEN
     341     CALL XERPRN (' *  ', -1, ' ', 72)
     342     CALL XERPRN (' ***', -1, 'END OF MESSAGE', 72)
     343     CALL XERPRN ('    ',  0, ' ', 72)
     344  ENDIF
     345  !
     346  !   IF THE ERROR IS NOT FATAL OR THE ERROR IS RECOVERABLE AND THE
     347  !   CONTROL FLAG IS SET FOR RECOVERY, THEN RETURN.
     348  !
     349   30   IF (LEVEL.LE.0 .OR. (LEVEL.EQ.1 .AND. MKNTRL.LE.1)) RETURN
     350  !
     351  !   THE PROGRAM WILL BE STOPPED DUE TO AN UNRECOVERED ERROR OR A
     352  !   FATAL ERROR.  PRINT THE REASON FOR THE ABORT AND THE ERROR
     353  !   SUMMARY IF THE CONTROL FLAG AND THE MAXIMUM ERROR COUNT PERMIT.
     354  !
     355  IF (LKNTRL.GT.0 .AND. KOUNT.LT.MAX(1,MAXMES)) THEN
     356     IF (LEVEL .EQ. 1) THEN
     357        CALL XERPRN &
     358              (' ***', -1, 'JOB ABORT DUE TO UNRECOVERED ERROR.', 72)
     359     ELSE
     360        CALL XERPRN(' ***', -1, 'JOB ABORT DUE TO FATAL ERROR.', 72)
     361     ENDIF
     362     CALL XERSVE (' ', ' ', ' ', -1, 0, 0, KDUMMY)
     363     CALL XERHLT (' ')
     364  ELSE
     365     CALL XERHLT (MESSG)
     366  ENDIF
     367  RETURN
     368END SUBROUTINE XERMSG

  • LMDZ6/trunk/libf/misc/xerprn.f90

    r5245 r5246  
    1 *DECK XERPRN
    2       SUBROUTINE XERPRN (PREFIX, NPREF, MESSG, NWRAP)
    3       IMPLICIT NONE
    4 C***BEGIN PROLOGUE  XERPRN
    5 C***SUBSIDIARY
    6 C***PURPOSE  Print error messages processed by XERMSG.
    7 C***LIBRARY   SLATEC (XERROR)
    8 C***CATEGORY  R3C
    9 C***TYPE      ALL (XERPRN-A)
    10 C***KEYWORDS  ERROR MESSAGES, PRINTING, XERROR
    11 C***AUTHOR  Fong, Kirby, (NMFECC at LLNL)
    12 C***DESCRIPTION
    13 C
    14 C This routine sends one or more lines to each of the (up to five)
    15 C logical units to which error messages are to be sent.  This routine
    16 C is called several times by XERMSG, sometimes with a single line to
    17 C print and sometimes with a (potentially very long) message that may
    18 C wrap around into multiple lines.
    19 C
    20 C PREFIX  Input argument of type CHARACTER.  This argument contains
    21 C         characters to be put at the beginning of each line before
    22 C         the body of the message.  No more than 16 characters of
    23 C         PREFIX will be used.
    24 C
    25 C NPREF   Input argument of type INTEGER.  This argument is the number
    26 C         of characters to use from PREFIX.  If it is negative, the
    27 C         intrinsic function LEN is used to determine its length.  If
    28 C         it is zero, PREFIX is not used.  If it exceeds 16 or if
    29 C         LEN(PREFIX) exceeds 16, only the first 16 characters will be
    30 C         used.  If NPREF is positive and the length of PREFIX is less
    31 C         than NPREF, a copy of PREFIX extended with blanks to length
    32 C         NPREF will be used.
    33 C
    34 C MESSG   Input argument of type CHARACTER.  This is the text of a
    35 C         message to be printed.  If it is a long message, it will be
    36 C         broken into pieces for printing on multiple lines.  Each line
    37 C         will start with the appropriate prefix and be followed by a
    38 C         piece of the message.  NWRAP is the number of characters per
    39 C         piece; that is, after each NWRAP characters, we break and
    40 C         start a new line.  In addition the characters '$$' embedded
    41 C         in MESSG are a sentinel for a new line.  The counting of
    42 C         characters up to NWRAP starts over for each new line.  The
    43 C         value of NWRAP typically used by XERMSG is 72 since many
    44 C         older error messages in the SLATEC Library are laid out to
    45 C         rely on wrap-around every 72 characters.
    46 C
    47 C NWRAP   Input argument of type INTEGER.  This gives the maximum size
    48 C         piece into which to break MESSG for printing on multiple
    49 C         lines.  An embedded '$$' ends a line, and the count restarts
    50 C         at the following character.  If a line break does not occur
    51 C         on a blank (it would split a word) that word is moved to the
    52 C         next line.  Values of NWRAP less than 16 will be treated as
    53 C         16.  Values of NWRAP greater than 132 will be treated as 132.
    54 C         The actual line length will be NPREF + NWRAP after NPREF has
    55 C         been adjusted to fall between 0 and 16 and NWRAP has been
    56 C         adjusted to fall between 16 and 132.
    57 C
    58 C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
    59 C                 Error-handling Package, SAND82-0800, Sandia
    60 C                 Laboratories, 1982.
    61 C***ROUTINES CALLED  I1MACH, XGETUA
    62 C***REVISION HISTORY  (YYMMDD)
    63 C   880621  DATE WRITTEN
    64 C   880708  REVISED AFTER THE SLATEC CML SUBCOMMITTEE MEETING OF
    65 C           JUNE 29 AND 30 TO CHANGE THE NAME TO XERPRN AND TO REWORK
    66 C           THE HANDLING OF THE NEW LINE SENTINEL TO BEHAVE LIKE THE
    67 C           SLASH CHARACTER IN FORMAT STATEMENTS.
    68 C   890706  REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO
    69 C           STREAMLINE THE CODING AND FIX A BUG THAT CAUSED EXTRA BLANK
    70 C           LINES TO BE PRINTED.
    71 C   890721  REVISED TO ADD A NEW FEATURE.  A NEGATIVE VALUE OF NPREF
    72 C           CAUSES LEN(PREFIX) TO BE USED AS THE LENGTH.
    73 C   891013  REVISED TO CORRECT ERROR IN CALCULATING PREFIX LENGTH.
    74 C   891214  Prologue converted to Version 4.0 format.  (WRB)
    75 C   900510  Added code to break messages between words.  (RWC)
    76 C   920501  Reformatted the REFERENCES section.  (WRB)
    77 C***END PROLOGUE  XERPRN
    78       CHARACTER*(*) PREFIX, MESSG
    79       INTEGER NPREF, NWRAP
    80       CHARACTER*148 CBUFF
    81       INTEGER IU(5), NUNIT
    82       CHARACTER*2 NEWLIN
    83       PARAMETER (NEWLIN = '$$')
    84       INTEGER N, I1MACH, I, LPREF, LWRAP, LENMSG, NEXTC
    85       INTEGER LPIECE, IDELTA
    86 C***FIRST EXECUTABLE STATEMENT  XERPRN
    87       CALL XGETUA(IU,NUNIT)
    88 C
    89 C       A ZERO VALUE FOR A LOGICAL UNIT NUMBER MEANS TO USE THE STANDARD
    90 C       ERROR MESSAGE UNIT INSTEAD.  I1MACH(4) RETRIEVES THE STANDARD
    91 C       ERROR MESSAGE UNIT.
    92 C
    93       N = I1MACH(4)
    94       DO 10 I=1,NUNIT
    95          IF (IU(I) .EQ. 0) IU(I) = N
    96    10 CONTINUE
    97 C
    98 C       LPREF IS THE LENGTH OF THE PREFIX.  THE PREFIX IS PLACED AT THE
    99 C       BEGINNING OF CBUFF, THE CHARACTER BUFFER, AND KEPT THERE DURING
    100 C       THE REST OF THIS ROUTINE.
    101 C
    102       IF ( NPREF .LT. 0 ) THEN
    103          LPREF = LEN(PREFIX)
    104       ELSE
    105          LPREF = NPREF
    106       ENDIF
    107       LPREF = MIN(16, LPREF)
    108       IF (LPREF .NE. 0) CBUFF(1:LPREF) = PREFIX
    109 C
    110 C       LWRAP IS THE MAXIMUM NUMBER OF CHARACTERS WE WANT TO TAKE AT ONE
    111 C       TIME FROM MESSG TO PRINT ON ONE LINE.
    112 C
    113       LWRAP = MAX(16, MIN(132, NWRAP))
    114 C
    115 C       SET LENMSG TO THE LENGTH OF MESSG, IGNORE ANY TRAILING BLANKS.
    116 C
    117       LENMSG = LEN(MESSG)
    118       N = LENMSG
    119       DO 20 I=1,N
    120          IF (MESSG(LENMSG:LENMSG) .NE. ' ') GO TO 30
    121          LENMSG = LENMSG - 1
    122    20 CONTINUE
    123    30 CONTINUE
    124 C
    125 C       IF THE MESSAGE IS ALL BLANKS, THEN PRINT ONE BLANK LINE.
    126 C
    127       IF (LENMSG .EQ. 0) THEN
    128          CBUFF(LPREF+1:LPREF+1) = ' '
    129          DO 40 I=1,NUNIT
    130             WRITE(IU(I), '(A)') CBUFF(1:LPREF+1)
    131    40    CONTINUE
    132          RETURN
    133       ENDIF
    134 C
    135 C       SET NEXTC TO THE POSITION IN MESSG WHERE THE NEXT SUBSTRING
    136 C       STARTS.  FROM THIS POSITION WE SCAN FOR THE NEW LINE SENTINEL.
    137 C       WHEN NEXTC EXCEEDS LENMSG, THERE IS NO MORE TO PRINT.
    138 C       WE LOOP BACK TO LABEL 50 UNTIL ALL PIECES HAVE BEEN PRINTED.
    139 C
    140 C       WE LOOK FOR THE NEXT OCCURRENCE OF THE NEW LINE SENTINEL.  THE
    141 C       INDEX INTRINSIC FUNCTION RETURNS ZERO IF THERE IS NO OCCURRENCE
    142 C       OR IF THE LENGTH OF THE FIRST ARGUMENT IS LESS THAN THE LENGTH
    143 C       OF THE SECOND ARGUMENT.
    144 C
    145 C       THERE ARE SEVERAL CASES WHICH SHOULD BE CHECKED FOR IN THE
    146 C       FOLLOWING ORDER.  WE ARE ATTEMPTING TO SET LPIECE TO THE NUMBER
    147 C       OF CHARACTERS THAT SHOULD BE TAKEN FROM MESSG STARTING AT
    148 C       POSITION NEXTC.
    149 C
    150 C       LPIECE .EQ. 0   THE NEW LINE SENTINEL DOES NOT OCCUR IN THE
    151 C                       REMAINDER OF THE CHARACTER STRING.  LPIECE
    152 C                       SHOULD BE SET TO LWRAP OR LENMSG+1-NEXTC,
    153 C                       WHICHEVER IS LESS.
    154 C
    155 C       LPIECE .EQ. 1   THE NEW LINE SENTINEL STARTS AT MESSG(NEXTC:
    156 C                       NEXTC).  LPIECE IS EFFECTIVELY ZERO, AND WE
    157 C                       PRINT NOTHING TO AVOID PRODUCING UNNECESSARY
    158 C                       BLANK LINES.  THIS TAKES CARE OF THE SITUATION
    159 C                       WHERE THE LIBRARY ROUTINE HAS A MESSAGE OF
    160 C                       EXACTLY 72 CHARACTERS FOLLOWED BY A NEW LINE
    161 C                       SENTINEL FOLLOWED BY MORE CHARACTERS.  NEXTC
    162 C                       SHOULD BE INCREMENTED BY 2.
    163 C
    164 C       LPIECE .GT. LWRAP+1  REDUCE LPIECE TO LWRAP.
    165 C
    166 C       ELSE            THIS LAST CASE MEANS 2 .LE. LPIECE .LE. LWRAP+1
    167 C                       RESET LPIECE = LPIECE-1.  NOTE THAT THIS
    168 C                       PROPERLY HANDLES THE END CASE WHERE LPIECE .EQ.
    169 C                       LWRAP+1.  THAT IS, THE SENTINEL FALLS EXACTLY
    170 C                       AT THE END OF A LINE.
    171 C
    172       NEXTC = 1
    173    50 LPIECE = INDEX(MESSG(NEXTC:LENMSG), NEWLIN)
    174       IF (LPIECE .EQ. 0) THEN
    175 C
    176 C       THERE WAS NO NEW LINE SENTINEL FOUND.
    177 C
    178          IDELTA = 0
    179          LPIECE = MIN(LWRAP, LENMSG+1-NEXTC)
    180          IF (LPIECE .LT. LENMSG+1-NEXTC) THEN
    181             DO 52 I=LPIECE+1,2,-1
    182                IF (MESSG(NEXTC+I-1:NEXTC+I-1) .EQ. ' ') THEN
    183                   LPIECE = I-1
    184                   IDELTA = 1
    185                   GOTO 54
    186                ENDIF
    187    52       CONTINUE
    188          ENDIF
    189    54    CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1)
    190          NEXTC = NEXTC + LPIECE + IDELTA
    191       ELSEIF (LPIECE .EQ. 1) THEN
    192 C
    193 C       WE HAVE A NEW LINE SENTINEL AT MESSG(NEXTC:NEXTC+1).
    194 C       DON'T PRINT A BLANK LINE.
    195 C
    196          NEXTC = NEXTC + 2
    197          GO TO 50
    198       ELSEIF (LPIECE .GT. LWRAP+1) THEN
    199 C
    200 C       LPIECE SHOULD BE SET DOWN TO LWRAP.
    201 C
    202          IDELTA = 0
    203          LPIECE = LWRAP
    204          DO 56 I=LPIECE+1,2,-1
    205             IF (MESSG(NEXTC+I-1:NEXTC+I-1) .EQ. ' ') THEN
    206                LPIECE = I-1
    207                IDELTA = 1
    208                GOTO 58
    209             ENDIF
    210    56    CONTINUE
    211    58    CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1)
    212          NEXTC = NEXTC + LPIECE + IDELTA
    213       ELSE
    214 C
    215 C       IF WE ARRIVE HERE, IT MEANS 2 .LE. LPIECE .LE. LWRAP+1.
    216 C       WE SHOULD DECREMENT LPIECE BY ONE.
    217 C
    218          LPIECE = LPIECE - 1
    219          CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1)
    220          NEXTC  = NEXTC + LPIECE + 2
    221       ENDIF
    222 C
    223 C       PRINT
    224 C
    225       DO 60 I=1,NUNIT
    226          WRITE(IU(I), '(A)') CBUFF(1:LPREF+LPIECE)
    227    60 CONTINUE
    228 C
    229       IF (NEXTC .LE. LENMSG) GO TO 50
    230       RETURN
    231       END
     1!DECK XERPRN
     2SUBROUTINE XERPRN (PREFIX, NPREF, MESSG, NWRAP)
     3  IMPLICIT NONE
     4  !***BEGIN PROLOGUE  XERPRN
     5  !***SUBSIDIARY
     6  !***PURPOSE  Print error messages processed by XERMSG.
     7  !***LIBRARY   SLATEC (XERROR)
     8  !***CATEGORY  R3C
     9  !***TYPE      ALL (XERPRN-A)
     10  !***KEYWORDS  ERROR MESSAGES, PRINTING, XERROR
     11  !***AUTHOR  Fong, Kirby, (NMFECC at LLNL)
     12  !***DESCRIPTION
     13  !
     14  ! This routine sends one or more lines to each of the (up to five)
     15  ! logical units to which error messages are to be sent.  This routine
     16  ! is called several times by XERMSG, sometimes with a single line to
     17  ! print and sometimes with a (potentially very long) message that may
     18  ! wrap around into multiple lines.
     19  !
     20  ! PREFIX  Input argument of type CHARACTER.  This argument contains
     21  !     characters to be put at the beginning of each line before
     22  !     the body of the message.  No more than 16 characters of
     23  !     PREFIX will be used.
     24  !
     25  ! NPREF   Input argument of type INTEGER.  This argument is the number
     26  !     of characters to use from PREFIX.  If it is negative, the
     27  !     intrinsic function LEN is used to determine its length.  If
     28  !         it is zero, PREFIX is not used.  If it exceeds 16 or if
     29  !     LEN(PREFIX) exceeds 16, only the first 16 characters will be
     30  !     used.  If NPREF is positive and the length of PREFIX is less
     31  !     than NPREF, a copy of PREFIX extended with blanks to length
     32  !     NPREF will be used.
     33  !
     34  ! MESSG   Input argument of type CHARACTER.  This is the text of a
     35  !     message to be printed.  If it is a long message, it will be
     36  !     broken into pieces for printing on multiple lines.  Each line
     37  !     will start with the appropriate prefix and be followed by a
     38  !     piece of the message.  NWRAP is the number of characters per
     39  !     piece; that is, after each NWRAP characters, we break and
     40  !     start a new line.  In addition the characters '$$' embedded
     41  !     in MESSG are a sentinel for a new line.  The counting of
     42  !     characters up to NWRAP starts over for each new line.  The
     43  !     value of NWRAP typically used by XERMSG is 72 since many
     44  !     older error messages in the SLATEC Library are laid out to
     45  !     rely on wrap-around every 72 characters.
     46  !
     47  ! NWRAP   Input argument of type INTEGER.  This gives the maximum size
     48  !     piece into which to break MESSG for printing on multiple
     49  !     lines.  An embedded '$$' ends a line, and the count restarts
     50  !     at the following character.  If a line break does not occur
     51  !     on a blank (it would split a word) that word is moved to the
     52  !     next line.  Values of NWRAP less than 16 will be treated as
     53  !     16.  Values of NWRAP greater than 132 will be treated as 132.
     54  !     The actual line length will be NPREF + NWRAP after NPREF has
     55  !     been adjusted to fall between 0 and 16 and NWRAP has been
     56  !     adjusted to fall between 16 and 132.
     57  !
     58  !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
     59  !             Error-handling Package, SAND82-0800, Sandia
     60  !             Laboratories, 1982.
     61  !***ROUTINES CALLED  I1MACH, XGETUA
     62  !***REVISION HISTORY  (YYMMDD)
     63  !   880621  DATE WRITTEN
     64  !   880708  REVISED AFTER THE SLATEC CML SUBCOMMITTEE MEETING OF
     65  !       JUNE 29 AND 30 TO CHANGE THE NAME TO XERPRN AND TO REWORK
     66  !       THE HANDLING OF THE NEW LINE SENTINEL TO BEHAVE LIKE THE
     67  !       SLASH CHARACTER IN FORMAT STATEMENTS.
     68  !   890706  REVISED WITH THE HELP OF FRED FRITSCH AND REG CLEMENS TO
     69  !       STREAMLINE THE CODING AND FIX A BUG THAT CAUSED EXTRA BLANK
     70  !       LINES TO BE PRINTED.
     71  !   890721  REVISED TO ADD A NEW FEATURE.  A NEGATIVE VALUE OF NPREF
     72  !       CAUSES LEN(PREFIX) TO BE USED AS THE LENGTH.
     73  !   891013  REVISED TO CORRECT ERROR IN CALCULATING PREFIX LENGTH.
     74  !   891214  Prologue converted to Version 4.0 format.  (WRB)
     75  !   900510  Added code to break messages between words.  (RWC)
     76  !   920501  Reformatted the REFERENCES section.  (WRB)
     77  !***END PROLOGUE  XERPRN
     78  CHARACTER(len=*) :: PREFIX, MESSG
     79  INTEGER :: NPREF, NWRAP
     80  CHARACTER(len=148) :: CBUFF
     81  INTEGER :: IU(5), NUNIT
     82  CHARACTER(len=2) :: NEWLIN
     83  PARAMETER (NEWLIN = '$$')
     84  INTEGER :: N, I1MACH, I, LPREF, LWRAP, LENMSG, NEXTC
     85  INTEGER :: LPIECE, IDELTA
     86  !***FIRST EXECUTABLE STATEMENT  XERPRN
     87  CALL XGETUA(IU,NUNIT)
     88  !
     89  !   A ZERO VALUE FOR A LOGICAL UNIT NUMBER MEANS TO USE THE STANDARD
     90  !   ERROR MESSAGE UNIT INSTEAD.  I1MACH(4) RETRIEVES THE STANDARD
     91  !   ERROR MESSAGE UNIT.
     92  !
     93  N = I1MACH(4)
     94  DO I=1,NUNIT
     95     IF (IU(I) .EQ. 0) IU(I) = N
     96  END DO
     97  !
     98  !   LPREF IS THE LENGTH OF THE PREFIX.  THE PREFIX IS PLACED AT THE
     99  !   BEGINNING OF CBUFF, THE CHARACTER BUFFER, AND KEPT THERE DURING
     100  !   THE REST OF THIS ROUTINE.
     101  !
     102  IF ( NPREF .LT. 0 ) THEN
     103     LPREF = LEN(PREFIX)
     104  ELSE
     105     LPREF = NPREF
     106  ENDIF
     107  LPREF = MIN(16, LPREF)
     108  IF (LPREF .NE. 0) CBUFF(1:LPREF) = PREFIX
     109  !
     110  !   LWRAP IS THE MAXIMUM NUMBER OF CHARACTERS WE WANT TO TAKE AT ONE
     111  !   TIME FROM MESSG TO PRINT ON ONE LINE.
     112  !
     113  LWRAP = MAX(16, MIN(132, NWRAP))
     114  !
     115  !   SET LENMSG TO THE LENGTH OF MESSG, IGNORE ANY TRAILING BLANKS.
     116  !
     117  LENMSG = LEN(MESSG)
     118  N = LENMSG
     119  DO I=1,N
     120     IF (MESSG(LENMSG:LENMSG) .NE. ' ') GO TO 30
     121     LENMSG = LENMSG - 1
     122  END DO
     123   30   CONTINUE
     124  !
     125  !   IF THE MESSAGE IS ALL BLANKS, THEN PRINT ONE BLANK LINE.
     126  !
     127  IF (LENMSG .EQ. 0) THEN
     128     CBUFF(LPREF+1:LPREF+1) = ' '
     129     DO I=1,NUNIT
     130        WRITE(IU(I), '(A)') CBUFF(1:LPREF+1)
     131     END DO
     132     RETURN
     133  ENDIF
     134  !
     135  !   SET NEXTC TO THE POSITION IN MESSG WHERE THE NEXT SUBSTRING
     136  !   STARTS.  FROM THIS POSITION WE SCAN FOR THE NEW LINE SENTINEL.
     137  !   WHEN NEXTC EXCEEDS LENMSG, THERE IS NO MORE TO PRINT.
     138  !   WE LOOP BACK TO LABEL 50 UNTIL ALL PIECES HAVE BEEN PRINTED.
     139  !
     140  !   WE LOOK FOR THE NEXT OCCURRENCE OF THE NEW LINE SENTINEL.  THE
     141  !   INDEX INTRINSIC FUNCTION RETURNS ZERO IF THERE IS NO OCCURRENCE
     142  !   OR IF THE LENGTH OF THE FIRST ARGUMENT IS LESS THAN THE LENGTH
     143  !   OF THE SECOND ARGUMENT.
     144  !
     145  !   THERE ARE SEVERAL CASES WHICH SHOULD BE CHECKED FOR IN THE
     146  !   FOLLOWING ORDER.  WE ARE ATTEMPTING TO SET LPIECE TO THE NUMBER
     147  !   OF CHARACTERS THAT SHOULD BE TAKEN FROM MESSG STARTING AT
     148  !   POSITION NEXTC.
     149  !
     150  !   LPIECE .EQ. 0   THE NEW LINE SENTINEL DOES NOT OCCUR IN THE
     151  !                   REMAINDER OF THE CHARACTER STRING.  LPIECE
     152  !                   SHOULD BE SET TO LWRAP OR LENMSG+1-NEXTC,
     153  !                   WHICHEVER IS LESS.
     154  !
     155  !   LPIECE .EQ. 1   THE NEW LINE SENTINEL STARTS AT MESSG(NEXTC:
     156  !                   NEXTC).  LPIECE IS EFFECTIVELY ZERO, AND WE
     157  !                   PRINT NOTHING TO AVOID PRODUCING UNNECESSARY
     158  !                   BLANK LINES.  THIS TAKES CARE OF THE SITUATION
     159  !                   WHERE THE LIBRARY ROUTINE HAS A MESSAGE OF
     160  !                   EXACTLY 72 CHARACTERS FOLLOWED BY A NEW LINE
     161  !                   SENTINEL FOLLOWED BY MORE CHARACTERS.  NEXTC
     162  !                   SHOULD BE INCREMENTED BY 2.
     163  !
     164  !   LPIECE .GT. LWRAP+1  REDUCE LPIECE TO LWRAP.
     165  !
     166  !   ELSE            THIS LAST CASE MEANS 2 .LE. LPIECE .LE. LWRAP+1
     167  !                   RESET LPIECE = LPIECE-1.  NOTE THAT THIS
     168  !                   PROPERLY HANDLES THE END CASE WHERE LPIECE .EQ.
     169  !                   LWRAP+1.  THAT IS, THE SENTINEL FALLS EXACTLY
     170  !                   AT THE END OF A LINE.
     171  !
     172  NEXTC = 1
     173   50   LPIECE = INDEX(MESSG(NEXTC:LENMSG), NEWLIN)
     174  IF (LPIECE .EQ. 0) THEN
     175  !
     176  !   THERE WAS NO NEW LINE SENTINEL FOUND.
     177  !
     178     IDELTA = 0
     179     LPIECE = MIN(LWRAP, LENMSG+1-NEXTC)
     180     IF (LPIECE .LT. LENMSG+1-NEXTC) THEN
     181        DO I=LPIECE+1,2,-1
     182           IF (MESSG(NEXTC+I-1:NEXTC+I-1) .EQ. ' ') THEN
     183              LPIECE = I-1
     184              IDELTA = 1
     185              GOTO 54
     186           ENDIF
     187        END DO
     188     ENDIF
     189   54   CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1)
     190     NEXTC = NEXTC + LPIECE + IDELTA
     191  ELSEIF (LPIECE .EQ. 1) THEN
     192  !
     193  !   WE HAVE A NEW LINE SENTINEL AT MESSG(NEXTC:NEXTC+1).
     194  !   DON'T PRINT A BLANK LINE.
     195  !
     196     NEXTC = NEXTC + 2
     197     GO TO 50
     198  ELSEIF (LPIECE .GT. LWRAP+1) THEN
     199  !
     200  !   LPIECE SHOULD BE SET DOWN TO LWRAP.
     201  !
     202     IDELTA = 0
     203     LPIECE = LWRAP
     204     DO I=LPIECE+1,2,-1
     205        IF (MESSG(NEXTC+I-1:NEXTC+I-1) .EQ. ' ') THEN
     206           LPIECE = I-1
     207           IDELTA = 1
     208           GOTO 58
     209        ENDIF
     210     END DO
     211   58   CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1)
     212     NEXTC = NEXTC + LPIECE + IDELTA
     213  ELSE
     214  !
     215  !   IF WE ARRIVE HERE, IT MEANS 2 .LE. LPIECE .LE. LWRAP+1.
     216  !   WE SHOULD DECREMENT LPIECE BY ONE.
     217  !
     218     LPIECE = LPIECE - 1
     219     CBUFF(LPREF+1:LPREF+LPIECE) = MESSG(NEXTC:NEXTC+LPIECE-1)
     220     NEXTC  = NEXTC + LPIECE + 2
     221  ENDIF
     222  !
     223  !   PRINT
     224  !
     225  DO I=1,NUNIT
     226     WRITE(IU(I), '(A)') CBUFF(1:LPREF+LPIECE)
     227  END DO
     228  !
     229  IF (NEXTC .LE. LENMSG) GO TO 50
     230  RETURN
     231END SUBROUTINE XERPRN

  • LMDZ6/trunk/libf/misc/xersve.f90

    r5245 r5246  
    1 *DECK XERSVE
    2       SUBROUTINE XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL,
    3      +   ICOUNT)
    4       IMPLICIT NONE
    5 C***BEGIN PROLOGUE  XERSVE
    6 C***SUBSIDIARY
    7 C***PURPOSE  Record that an error has occurred.
    8 C***LIBRARY   SLATEC (XERROR)
    9 C***CATEGORY  R3
    10 C***TYPE      ALL (XERSVE-A)
    11 C***KEYWORDS  ERROR, XERROR
    12 C***AUTHOR  Jones, R. E., (SNLA)
    13 C***DESCRIPTION
    14 C
    15 C *Usage:
    16 C
    17 C        INTEGER  KFLAG, NERR, LEVEL, ICOUNT
    18 C        CHARACTER * (len) LIBRAR, SUBROU, MESSG
    19 C
    20 C        CALL XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL, ICOUNT)
    21 C
    22 C *Arguments:
    23 C
    24 C        LIBRAR :IN    is the library that the message is from.
    25 C        SUBROU :IN    is the subroutine that the message is from.
    26 C        MESSG  :IN    is the message to be saved.
    27 C        KFLAG  :IN    indicates the action to be performed.
    28 C                      when KFLAG > 0, the message in MESSG is saved.
    29 C                      when KFLAG=0 the tables will be dumped and
    30 C                      cleared.
    31 C                      when KFLAG < 0, the tables will be dumped and
    32 C                      not cleared.
    33 C        NERR   :IN    is the error number.
    34 C        LEVEL  :IN    is the error severity.
    35 C        ICOUNT :OUT   the number of times this message has been seen,
    36 C                      or zero if the table has overflowed and does not
    37 C                      contain this message specifically.  When KFLAG=0,
    38 C                      ICOUNT will not be altered.
    39 C
    40 C *Description:
    41 C
    42 C   Record that this error occurred and possibly dump and clear the
    43 C   tables.
    44 C
    45 C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
    46 C                 Error-handling Package, SAND82-0800, Sandia
    47 C                 Laboratories, 1982.
    48 C***ROUTINES CALLED  I1MACH, XGETUA
    49 C***REVISION HISTORY  (YYMMDD)
    50 C   800319  DATE WRITTEN
    51 C   861211  REVISION DATE from Version 3.2
    52 C   891214  Prologue converted to Version 4.0 format.  (BAB)
    53 C   900413  Routine modified to remove reference to KFLAG.  (WRB)
    54 C   900510  Changed to add LIBRARY NAME and SUBROUTINE to calling
    55 C           sequence, use IF-THEN-ELSE, make number of saved entries
    56 C           easily changeable, changed routine name from XERSAV to
    57 C           XERSVE.  (RWC)
    58 C   910626  Added LIBTAB and SUBTAB to SAVE statement.  (BKS)
    59 C   920501  Reformatted the REFERENCES section.  (WRB)
    60 C***END PROLOGUE  XERSVE
    61       INTEGER,PARAMETER :: LENTAB=10
    62       INTEGER LUN(5)
    63       CHARACTER*(*) LIBRAR, SUBROU, MESSG
    64       CHARACTER*8 LIBTAB(LENTAB), SUBTAB(LENTAB), LIB, SUB
    65       CHARACTER*20 MESTAB(LENTAB), MES
    66       DIMENSION NERTAB(LENTAB), LEVTAB(LENTAB), KOUNT(LENTAB)
    67       SAVE LIBTAB, SUBTAB, MESTAB, NERTAB, LEVTAB, KOUNT, KOUNTX, NMSG
    68       DATA KOUNTX/0/, NMSG/0/
    69       INTEGER NERR,LEVEL,KONTRL
    70       INTEGER NERTAB, LEVTAB, KOUNT, KOUNTX, NMSG
    71       INTEGER KFLAG, ICOUNT, NUNIT, KUNIT, IUNIT, I1MACH, I
    72 C***FIRST EXECUTABLE STATEMENT  XERSVE
    73 C
    74       IF (KFLAG.LE.0) THEN
    75 C
    76 C        Dump the table.
    77 C
    78          IF (NMSG.EQ.0) RETURN
    79 C
    80 C        Print to each unit.
    81 C
    82          CALL XGETUA (LUN, NUNIT)
    83          DO 20 KUNIT = 1,NUNIT
    84             IUNIT = LUN(KUNIT)
    85             IF (IUNIT.EQ.0) IUNIT = I1MACH(4)
    86 C
    87 C           Print the table header.
    88 C
    89             WRITE (IUNIT,9000)
    90 C
    91 C           Print body of table.
    92 C
    93             DO 10 I = 1,NMSG
    94                WRITE (IUNIT,9010) LIBTAB(I), SUBTAB(I), MESTAB(I),
    95      *            NERTAB(I),LEVTAB(I),KOUNT(I)
    96    10       CONTINUE
    97 C
    98 C           Print number of other errors.
    99 C
    100             IF (KOUNTX.NE.0) WRITE (IUNIT,9020) KOUNTX
    101             WRITE (IUNIT,9030)
    102    20    CONTINUE
    103 C
    104 C        Clear the error tables.
    105 C
    106          IF (KFLAG.EQ.0) THEN
    107             NMSG = 0
    108             KOUNTX = 0
    109          ENDIF
    110       ELSE
    111 C
    112 C        PROCESS A MESSAGE...
    113 C        SEARCH FOR THIS MESSG, OR ELSE AN EMPTY SLOT FOR THIS MESSG,
    114 C        OR ELSE DETERMINE THAT THE ERROR TABLE IS FULL.
    115 C
    116          LIB = LIBRAR
    117          SUB = SUBROU
    118          MES = MESSG
    119          DO 30 I = 1,NMSG
    120             IF (LIB.EQ.LIBTAB(I) .AND. SUB.EQ.SUBTAB(I) .AND.
    121      *         MES.EQ.MESTAB(I) .AND. NERR.EQ.NERTAB(I) .AND.
    122      *         LEVEL.EQ.LEVTAB(I)) THEN
    123                   KOUNT(I) = KOUNT(I) + 1
    124                   ICOUNT = KOUNT(I)
    125                   RETURN
    126             ENDIF
    127    30    CONTINUE
    128 C
    129          IF (NMSG.LT.LENTAB) THEN
    130 C
    131 C           Empty slot found for new message.
    132 C
    133             NMSG = NMSG + 1
    134             LIBTAB(I) = LIB
    135             SUBTAB(I) = SUB
    136             MESTAB(I) = MES
    137             NERTAB(I) = NERR
    138             LEVTAB(I) = LEVEL
    139             KOUNT (I) = 1
    140             ICOUNT    = 1
    141          ELSE
    142 C
    143 C           Table is full.
    144 C
    145             KOUNTX = KOUNTX+1
    146             ICOUNT = 0
    147          ENDIF
    148       ENDIF
    149       RETURN
    150 C
    151 C    Formats.
    152 C
    153  9000 FORMAT ('0          ERROR MESSAGE SUMMARY' /
    154      +   ' LIBRARY    SUBROUTINE MESSAGE START             NERR',
    155      +   '     LEVEL     COUNT')
    156  9010 FORMAT (1X,A,3X,A,3X,A,3I10)
    157  9020 FORMAT ('0OTHER ERRORS NOT INDIVIDUALLY TABULATED = ', I10)
    158  9030 FORMAT (1X)
    159       END
     1!DECK XERSVE
     2SUBROUTINE XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL, &
     3        ICOUNT)
     4  IMPLICIT NONE
     5  !***BEGIN PROLOGUE  XERSVE
     6  !***SUBSIDIARY
     7  !***PURPOSE  Record that an error has occurred.
     8  !***LIBRARY   SLATEC (XERROR)
     9  !***CATEGORY  R3
     10  !***TYPE      ALL (XERSVE-A)
     11  !***KEYWORDS  ERROR, XERROR
     12  !***AUTHOR  Jones, R. E., (SNLA)
     13  !***DESCRIPTION
     14  !
     15  ! *Usage:
     16  !
     17  !    INTEGER  KFLAG, NERR, LEVEL, ICOUNT
     18  !    CHARACTER * (len) LIBRAR, SUBROU, MESSG
     19  !
     20  !    CALL XERSVE (LIBRAR, SUBROU, MESSG, KFLAG, NERR, LEVEL, ICOUNT)
     21  !
     22  ! *Arguments:
     23  !
     24  !    LIBRAR :IN    is the library that the message is from.
     25  !    SUBROU :IN    is the subroutine that the message is from.
     26  !    MESSG  :IN    is the message to be saved.
     27  !    KFLAG  :IN    indicates the action to be performed.
     28  !                  when KFLAG > 0, the message in MESSG is saved.
     29  !                  when KFLAG=0 the tables will be dumped and
     30  !                  cleared.
     31  !                  when KFLAG < 0, the tables will be dumped and
     32  !                  not cleared.
     33  !    NERR   :IN    is the error number.
     34  !    LEVEL  :IN    is the error severity.
     35  !    ICOUNT :OUT   the number of times this message has been seen,
     36  !                  or zero if the table has overflowed and does not
     37  !                  contain this message specifically.  When KFLAG=0,
     38  !                  ICOUNT will not be altered.
     39  !
     40  ! *Description:
     41  !
     42  !   Record that this error occurred and possibly dump and clear the
     43  !   tables.
     44  !
     45  !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
     46  !             Error-handling Package, SAND82-0800, Sandia
     47  !             Laboratories, 1982.
     48  !***ROUTINES CALLED  I1MACH, XGETUA
     49  !***REVISION HISTORY  (YYMMDD)
     50  !   800319  DATE WRITTEN
     51  !   861211  REVISION DATE from Version 3.2
     52  !   891214  Prologue converted to Version 4.0 format.  (BAB)
     53  !   900413  Routine modified to remove reference to KFLAG.  (WRB)
     54  !   900510  Changed to add LIBRARY NAME and SUBROUTINE to calling
     55  !       sequence, use IF-THEN-ELSE, make number of saved entries
     56  !       easily changeable, changed routine name from XERSAV to
     57  !       XERSVE.  (RWC)
     58  !   910626  Added LIBTAB and SUBTAB to SAVE statement.  (BKS)
     59  !   920501  Reformatted the REFERENCES section.  (WRB)
     60  !***END PROLOGUE  XERSVE
     61  INTEGER,PARAMETER :: LENTAB=10
     62  INTEGER :: LUN(5)
     63  CHARACTER(len=*) :: LIBRAR, SUBROU, MESSG
     64  CHARACTER(len=8) :: LIBTAB(LENTAB), SUBTAB(LENTAB), LIB, SUB
     65  CHARACTER(len=20) :: MESTAB(LENTAB), MES
     66  DIMENSION NERTAB(LENTAB), LEVTAB(LENTAB), KOUNT(LENTAB)
     67  SAVE LIBTAB, SUBTAB, MESTAB, NERTAB, LEVTAB, KOUNT, KOUNTX, NMSG
     68  DATA KOUNTX/0/, NMSG/0/
     69  INTEGER :: NERR,LEVEL,KONTRL
     70  INTEGER :: NERTAB, LEVTAB, KOUNT, KOUNTX, NMSG
     71  INTEGER :: KFLAG, ICOUNT, NUNIT, KUNIT, IUNIT, I1MACH, I
     72  !***FIRST EXECUTABLE STATEMENT  XERSVE
     73  !
     74  IF (KFLAG.LE.0) THEN
     75  !
     76  !    Dump the table.
     77  !
     78     IF (NMSG.EQ.0) RETURN
     79  !
     80  !    Print to each unit.
     81  !
     82     CALL XGETUA (LUN, NUNIT)
     83     DO KUNIT = 1,NUNIT
     84        IUNIT = LUN(KUNIT)
     85        IF (IUNIT.EQ.0) IUNIT = I1MACH(4)
     86  !
     87  !       Print the table header.
     88  !
     89        WRITE (IUNIT,9000)
     90  !
     91  !       Print body of table.
     92  !
     93        DO I = 1,NMSG
     94           WRITE (IUNIT,9010) LIBTAB(I), SUBTAB(I), MESTAB(I), &
     95                 NERTAB(I),LEVTAB(I),KOUNT(I)
     96        END DO
     97  !
     98  !       Print number of other errors.
     99  !
     100        IF (KOUNTX.NE.0) WRITE (IUNIT,9020) KOUNTX
     101        WRITE (IUNIT,9030)
     102     END DO
     103  !
     104  !    Clear the error tables.
     105  !
     106     IF (KFLAG.EQ.0) THEN
     107        NMSG = 0
     108        KOUNTX = 0
     109     ENDIF
     110  ELSE
     111  !
     112  !    PROCESS A MESSAGE...
     113  !    SEARCH FOR THIS MESSG, OR ELSE AN EMPTY SLOT FOR THIS MESSG,
     114  !    OR ELSE DETERMINE THAT THE ERROR TABLE IS FULL.
     115  !
     116     LIB = LIBRAR
     117     SUB = SUBROU
     118     MES = MESSG
     119     DO I = 1,NMSG
     120        IF (LIB.EQ.LIBTAB(I) .AND. SUB.EQ.SUBTAB(I) .AND. &
     121              MES.EQ.MESTAB(I) .AND. NERR.EQ.NERTAB(I) .AND. &
     122              LEVEL.EQ.LEVTAB(I)) THEN
     123              KOUNT(I) = KOUNT(I) + 1
     124              ICOUNT = KOUNT(I)
     125              RETURN
     126        ENDIF
     127     END DO
     128  !
     129     IF (NMSG.LT.LENTAB) THEN
     130  !
     131  !       Empty slot found for new message.
     132  !
     133        NMSG = NMSG + 1
     134        LIBTAB(I) = LIB
     135        SUBTAB(I) = SUB
     136        MESTAB(I) = MES
     137        NERTAB(I) = NERR
     138        LEVTAB(I) = LEVEL
     139        KOUNT (I) = 1
     140        ICOUNT    = 1
     141     ELSE
     142  !
     143  !       Table is full.
     144  !
     145        KOUNTX = KOUNTX+1
     146        ICOUNT = 0
     147     ENDIF
     148  ENDIF
     149  RETURN
     150  !
     151  ! Formats.
     152  !
     153 9000   FORMAT ('0          ERROR MESSAGE SUMMARY' / &
     154              ' LIBRARY    SUBROUTINE MESSAGE START             NERR', &
     155              '     LEVEL     COUNT')
     156 9010   FORMAT (1X,A,3X,A,3X,A,3I10)
     157 9020   FORMAT ('0OTHER ERRORS NOT INDIVIDUALLY TABULATED = ', I10)
     158 9030   FORMAT (1X)
     159END SUBROUTINE XERSVE

  • LMDZ6/trunk/libf/misc/xgetua.f90

    r5245 r5246  
    1 *DECK XGETUA
    2       SUBROUTINE XGETUA (IUNITA, N)
    3       IMPLICIT NONE
    4 C***BEGIN PROLOGUE  XGETUA
    5 C***PURPOSE  Return unit number(s) to which error messages are being
    6 C            sent.
    7 C***LIBRARY   SLATEC (XERROR)
    8 C***CATEGORY  R3C
    9 C***TYPE      ALL (XGETUA-A)
    10 C***KEYWORDS  ERROR, XERROR
    11 C***AUTHOR  Jones, R. E., (SNLA)
    12 C***DESCRIPTION
    13 C
    14 C    Abstract
    15 C        XGETUA may be called to determine the unit number or numbers
    16 C        to which error messages are being sent.
    17 C        These unit numbers may have been set by a call to XSETUN,
    18 C        or a call to XSETUA, or may be a default value.
    19 C
    20 C    Description of Parameters
    21 C      --Output--
    22 C        IUNIT - an array of one to five unit numbers, depending
    23 C                on the value of N.  A value of zero refers to the
    24 C                default unit, as defined by the I1MACH machine
    25 C                constant routine.  Only IUNIT(1),...,IUNIT(N) are
    26 C                defined by XGETUA.  The values of IUNIT(N+1),...,
    27 C                IUNIT(5) are not defined (for N .LT. 5) or altered
    28 C                in any way by XGETUA.
    29 C        N     - the number of units to which copies of the
    30 C                error messages are being sent.  N will be in the
    31 C                range from 1 to 5.
    32 C
    33 C***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
    34 C                 Error-handling Package, SAND82-0800, Sandia
    35 C                 Laboratories, 1982.
    36 C***ROUTINES CALLED  J4SAVE
    37 C***REVISION HISTORY  (YYMMDD)
    38 C   790801  DATE WRITTEN
    39 C   861211  REVISION DATE from Version 3.2
    40 C   891214  Prologue converted to Version 4.0 format.  (BAB)
    41 C   920501  Reformatted the REFERENCES section.  (WRB)
    42 C***END PROLOGUE  XGETUA
    43       DIMENSION IUNITA(5)
    44       INTEGER IUNITA, N, J4SAVE, INDEX, I
    45 C***FIRST EXECUTABLE STATEMENT  XGETUA
    46       N = J4SAVE(5,0,.FALSE.)
    47       DO 30 I=1,N
    48          INDEX = I+4
    49          IF (I.EQ.1) INDEX = 3
    50          IUNITA(I) = J4SAVE(INDEX,0,.FALSE.)
    51    30 CONTINUE
    52       RETURN
    53       END
     1!DECK XGETUA
     2SUBROUTINE XGETUA (IUNITA, N)
     3  IMPLICIT NONE
     4  !***BEGIN PROLOGUE  XGETUA
     5  !***PURPOSE  Return unit number(s) to which error messages are being
     6         ! sent.
     7  !***LIBRARY   SLATEC (XERROR)
     8  !***CATEGORY  R3C
     9  !***TYPE      ALL (XGETUA-A)
     10  !***KEYWORDS  ERROR, XERROR
     11  !***AUTHOR  Jones, R. E., (SNLA)
     12  !***DESCRIPTION
     13  !
     14  ! Abstract
     15  !    XGETUA may be called to determine the unit number or numbers
     16  !    to which error messages are being sent.
     17  !    These unit numbers may have been set by a call to XSETUN,
     18  !    or a call to XSETUA, or may be a default value.
     19  !
     20  ! Description of Parameters
     21  !  --Output--
     22  !    IUNIT - an array of one to five unit numbers, depending
     23  !            on the value of N.  A value of zero refers to the
     24  !            default unit, as defined by the I1MACH machine
     25  !            constant routine.  Only IUNIT(1),...,IUNIT(N) are
     26  !            defined by XGETUA.  The values of IUNIT(N+1),...,
     27  !            IUNIT(5) are not defined (for N .LT. 5) or altered
     28  !            in any way by XGETUA.
     29  !    N     - the number of units to which copies of the
     30  !            error messages are being sent.  N will be in the
     31  !            range from 1 to 5.
     32  !
     33  !***REFERENCES  R. E. Jones and D. K. Kahaner, XERROR, the SLATEC
     34  !             Error-handling Package, SAND82-0800, Sandia
     35  !             Laboratories, 1982.
     36  !***ROUTINES CALLED  J4SAVE
     37  !***REVISION HISTORY  (YYMMDD)
     38  !   790801  DATE WRITTEN
     39  !   861211  REVISION DATE from Version 3.2
     40  !   891214  Prologue converted to Version 4.0 format.  (BAB)
     41  !   920501  Reformatted the REFERENCES section.  (WRB)
     42  !***END PROLOGUE  XGETUA
     43  DIMENSION IUNITA(5)
     44  INTEGER :: IUNITA, N, J4SAVE, INDEX, I
     45  !***FIRST EXECUTABLE STATEMENT  XGETUA
     46  N = J4SAVE(5,0,.FALSE.)
     47  DO I=1,N
     48     INDEX = I+4
     49     IF (I.EQ.1) INDEX = 3
     50     IUNITA(I) = J4SAVE(INDEX,0,.FALSE.)
     51  END DO
     52  RETURN
     53END SUBROUTINE XGETUA
Note: See TracChangeset for help on using the changeset viewer.