Ignore:
Timestamp:
Jul 23, 2024, 7:14:34 PM (12 months ago)
Author:
abarral
Message:

Replace 1DUTILS.h by module lmdz_1dutils.f90
Replace 1DConv.h by module lmdz_old_1dconv.f90 (it's only used by old_* files)
Convert *.F to *.f90
Fix gradsdef.h formatting
Remove unnecessary "RETURN" at the end of functions/subroutines

Location:
LMDZ6/branches/Amaury_dev/libf/misc
Files:
23 moved

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Amaury_dev/libf/misc/cbrt.f90

    r5104 r5105  
    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
    119
     10END FUNCTION cbrt
     11

  • LMDZ6/branches/Amaury_dev/libf/misc/chfev.f90

    r5104 r5105  
    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 < 1)  GO TO 5001
    103       H = X2 - X1
    104       IF (H == 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 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<XMI )  NEXT(1) = NEXT(1) + 1
    131          IF ( X>XMA )  NEXT(2) = NEXT(2) + 1
    132 C        (NOTE REDUNDANCY--IF EITHER CONDITION IS TRUE, OTHER IS FALSE.)
    133       END DO
    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 < 1)  GO TO 5001
     103  H = X2 - X1
     104  IF (H == 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<XMI )  NEXT(1) = NEXT(1) + 1
     131     IF ( X>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/branches/Amaury_dev/libf/misc/fdump.f90

    r5104 r5105  
    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
     31END SUBROUTINE FDUMP

  • LMDZ6/branches/Amaury_dev/libf/misc/formcoord.f90

    r5104 r5105  
    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<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<dxmin) dxmin=dx
    35          enddo
     27  if (n<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<dxmin) dxmin=dx
     35     enddo
    3636
    37          ndec=-log10(dxmin)+2
    38          if(mod(n,6)==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)==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')
    5352
    54       end
     53
     54end subroutine formcoord

  • LMDZ6/branches/Amaury_dev/libf/misc/i1mach.f90

    r5104 r5105  
    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 < 1  .OR.  I > 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 < 1  .OR.  I > 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/branches/Amaury_dev/libf/misc/ismax.f90

    r5104 r5105  
    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 i=1,n-1
    15        ix=ix+incx
    16        if(sx(ix)>sxmax) then
    17          sxmax=sx(ix)
    18          ismax=i+1
    19        endif
    20       END DO
    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)>sxmax) then
     17     sxmax=sx(ix)
     18     ismax=i+1
     19   endif
     20  END DO
     21  !
    2422
     23end function ismax
     24

  • LMDZ6/branches/Amaury_dev/libf/misc/ismin.f90

    r5104 r5105  
    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)<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)<sxmin) then
     17         sxmin=sx(ix)
     18         ismin=i+1
     19     endif
     20  ENDDO
     21  !
     22
     23end function ismin
     24!

  • LMDZ6/branches/Amaury_dev/libf/misc/j4save.f90

    r5104 r5105  
    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
     67END FUNCTION J4SAVE

  • LMDZ6/branches/Amaury_dev/libf/misc/juldate.f90

    r5104 r5105  
    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 <= 2) then
    21             year=year-1.
    22             rmon=rmon+12.
    23         endif
    24         cf=year+(rmon/100.)+(ojou/10000.)
    25         if (cf >= 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 <= 2) then
     21        year=year-1.
     22        rmon=rmon+12.
     23    endif
     24    cf=year+(rmon/100.)+(ojou/10000.)
     25    if (cf >= 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
     37end subroutine juldate
    3838
    3939

  • LMDZ6/branches/Amaury_dev/libf/misc/minmax.f90

    r5104 r5105  
    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
    2321
     22END SUBROUTINE minmax
     23

  • LMDZ6/branches/Amaury_dev/libf/misc/minmax2.f90

    r5104 r5105  
    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
     20END SUBROUTINE minmax2

  • LMDZ6/branches/Amaury_dev/libf/misc/pchdf.f90

    r5104 r5105  
    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 < 3)  GO TO 5001
    75 C
    76 C  COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL.
    77 C
    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 C
    84 C  EVALUATE DERIVATIVE AT X(K).
    85 C
    86       VALUE = S(1)
    87       DO I = 2, K-1
    88          VALUE = S(I) + VALUE*(X(K)-X(I))
    89       END DO
    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 < 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/branches/Amaury_dev/libf/misc/pchfe.f90

    r5104 r5105  
    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<2 )  GO TO 5001
    146       IF ( INCFD<1 )  GO TO 5002
    147       DO I = 2, N
    148          IF ( X(I)<=X(I-1) )  GO TO 5003
    149       END DO
    150 C
    151 C  FUNCTION DEFINITION IS OK, GO ON.
    152 C
    153     5 CONTINUE
    154       IF ( NE<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 > NE)  GO TO 5000
    167 C
    168 C    LOCATE ALL POINTS IN INTERVAL.
    169 C
    170          DO J = JFIRST, NE
    171             IF (XE(J) >= X(IR))  GO TO 30
    172       END DO
    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 == 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 == 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 < 0)  GO TO 5005
    195 C
    196          IF (NEXT(2) == 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 < 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) == 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 > 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 I = JFIRST, J-1
    231                   IF (XE(I) < X(IR-1))  GO TO 45
    232       END DO
    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 I = 1, IR-1
    243                   IF (XE(J) < X(I)) GO TO 47
    244       END DO
    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 <= 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<2 )  GO TO 5001
     146  IF ( INCFD<1 )  GO TO 5002
     147  DO I = 2, N
     148     IF ( X(I)<=X(I-1) )  GO TO 5003
     149  END DO
     150  !
     151  !  FUNCTION DEFINITION IS OK, GO ON.
     152  !
     153    5   CONTINUE
     154  IF ( NE<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 > NE)  GO TO 5000
     167  !
     168  ! LOCATE ALL POINTS IN INTERVAL.
     169  !
     170     DO J = JFIRST, NE
     171        IF (XE(J) >= 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 == 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 == 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 < 0)  GO TO 5005
     195  !
     196     IF (NEXT(2) == 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 < 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) == 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 > 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) < 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) < 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 <= 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/branches/Amaury_dev/libf/misc/pchsp.f90

    r5104 r5105  
    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<2 )  GO TO 5001
    165       IF ( INCFD<1 )  GO TO 5002
    166       DO J = 2, N
    167          IF ( X(J)<=X(J-1) )  GO TO 5003
    168       END DO
    169 C
    170       IBEG = IC(1)
    171       IEND = IC(2)
    172       IERR = 0
    173       IF ( (IBEG<0).OR.(IBEG>4) )  IERR = IERR - 1
    174       IF ( (IEND<0).OR.(IEND>4) )  IERR = IERR - 2
    175       IF ( IERR<0 )  GO TO 5004
    176 C
    177 C  FUNCTION DEFINITION IS OK -- GO ON.
    178 C
    179       IF ( NWK < 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 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 C
    188 C  SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL.
    189 C
    190       IF ( IBEG>N )  IBEG = 0
    191       IF ( IEND>N )  IEND = 0
    192 C
    193 C  SET UP FOR BOUNDARY CONDITIONS.
    194 C
    195       IF ( (IBEG==1).OR.(IBEG==2) )  THEN
    196          D(1,1) = VC(1)
    197       ELSE IF (IBEG > 2)  THEN
    198 C        PICK UP FIRST IBEG POINTS, IN REVERSE ORDER.
    199          DO 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 < IBEG)  STEMP(J) = WK(2,INDEX)
    204       END DO
    205 C                --------------------------------
    206          D(1,1) = PCHDF (IBEG, XTEMP, STEMP, IERR)
    207 C                --------------------------------
    208          IF (IERR /= 0)  GO TO 5009
    209          IBEG = 1
    210       ENDIF
    211 C
    212       IF ( (IEND==1).OR.(IEND==2) )  THEN
    213          D(1,N) = VC(2)
    214       ELSE IF (IEND > 2)  THEN
    215 C        PICK UP LAST IEND POINTS.
    216          DO 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 < IEND)  STEMP(J) = WK(2,INDEX+1)
    221       END DO
    222 C                --------------------------------
    223          D(1,N) = PCHDF (IEND, XTEMP, STEMP, IERR)
    224 C                --------------------------------
    225          IF (IERR /= 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 == 0)  THEN
    240          IF (N == 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 == 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 > 1)  THEN
    269          DO J=2,NM1
    270             IF (WK(2,J-1) == 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 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 == 1)  GO TO 30
    285 C
    286       IF (IEND == 0)  THEN
    287          IF (N==2 .AND. IBEG==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==2) .OR. (N==3 .AND. IBEG==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) == 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) == 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) == 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) == 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 J=NM1,1,-1
    327          IF (WK(2,J) == 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 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<2 )  GO TO 5001
     165  IF ( INCFD<1 )  GO TO 5002
     166  DO J = 2, N
     167     IF ( X(J)<=X(J-1) )  GO TO 5003
     168  END DO
     169  !
     170  IBEG = IC(1)
     171  IEND = IC(2)
     172  IERR = 0
     173  IF ( (IBEG<0).OR.(IBEG>4) )  IERR = IERR - 1
     174  IF ( (IEND<0).OR.(IEND>4) )  IERR = IERR - 2
     175  IF ( IERR<0 )  GO TO 5004
     176  !
     177  !  FUNCTION DEFINITION IS OK -- GO ON.
     178  !
     179  IF ( NWK < 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>N )  IBEG = 0
     191  IF ( IEND>N )  IEND = 0
     192  !
     193  !  SET UP FOR BOUNDARY CONDITIONS.
     194  !
     195  IF ( (IBEG==1).OR.(IBEG==2) )  THEN
     196     D(1,1) = VC(1)
     197  ELSE IF (IBEG > 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 < IBEG)  STEMP(J) = WK(2,INDEX)
     204  END DO
     205              ! --------------------------------
     206     D(1,1) = PCHDF (IBEG, XTEMP, STEMP, IERR)
     207              ! --------------------------------
     208     IF (IERR /= 0)  GO TO 5009
     209     IBEG = 1
     210  ENDIF
     211  !
     212  IF ( (IEND==1).OR.(IEND==2) )  THEN
     213     D(1,N) = VC(2)
     214  ELSE IF (IEND > 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 < IEND)  STEMP(J) = WK(2,INDEX+1)
     221  END DO
     222              ! --------------------------------
     223     D(1,N) = PCHDF (IEND, XTEMP, STEMP, IERR)
     224              ! --------------------------------
     225     IF (IERR /= 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 == 0)  THEN
     240     IF (N == 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 == 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 > 1)  THEN
     269     DO J=2,NM1
     270        IF (WK(2,J-1) == 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 == 1)  GO TO 30
     285  !
     286  IF (IEND == 0)  THEN
     287     IF (N==2 .AND. IBEG==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==2) .OR. (N==3 .AND. IBEG==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) == 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) == 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) == 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) == 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) == 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/branches/Amaury_dev/libf/misc/q_sat.f90

    r5104 r5105  
    44
    55
    6       SUBROUTINE q_sat(np,temp,pres,qsat)
     6SUBROUTINE q_sat(np,temp,pres,qsat)
    77
    8       IMPLICIT none
    9 !======================================================================
    10 ! Autheur(s): Z.X. Li (LMD/CNRS)
    11 !  reecriture vectorisee par F. Hourdin.
    12 ! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
    13 !======================================================================
    14 ! Arguments:
    15 ! kelvin---input-R: temperature en Kelvin
    16 ! millibar--input-R: pression en mb
     8  IMPLICIT none
     9  !======================================================================
     10  ! Autheur(s): Z.X. Li (LMD/CNRS)
     11  !  reecriture vectorisee par F. Hourdin.
     12  ! Objet: calculer la vapeur d'eau saturante (formule Centre Euro.)
     13  !======================================================================
     14  ! Arguments:
     15  ! kelvin---input-R: temperature en Kelvin
     16  ! millibar--input-R: pression en mb
    1717
    18 ! q_sat----output-R: vapeur d'eau saturante en kg/kg
    19 !======================================================================
     18  ! q_sat----output-R: vapeur d'eau saturante en kg/kg
     19  !======================================================================
    2020
    21       integer np
    22       REAL temp(np),pres(np),qsat(np)
     21  integer :: np
     22  REAL :: temp(np),pres(np),qsat(np)
    2323
    24       REAL r2es
    25       PARAMETER (r2es=611.14 *18.0153/28.9644)
     24  REAL :: r2es
     25  PARAMETER (r2es=611.14 *18.0153/28.9644)
    2626
    27       REAL r3les, r3ies, r3es
    28       PARAMETER (R3LES=17.269)
    29       PARAMETER (R3IES=21.875)
     27  REAL :: r3les, r3ies, r3es
     28  PARAMETER (R3LES=17.269)
     29  PARAMETER (R3IES=21.875)
    3030
    31       REAL r4les, r4ies, r4es
    32       PARAMETER (R4LES=35.86)
    33       PARAMETER (R4IES=7.66)
     31  REAL :: r4les, r4ies, r4es
     32  PARAMETER (R4LES=35.86)
     33  PARAMETER (R4IES=7.66)
    3434
    35       REAL rtt
    36       PARAMETER (rtt=273.16)
     35  REAL :: rtt
     36  PARAMETER (rtt=273.16)
    3737
    38       REAL retv
    39       PARAMETER (retv=28.9644/18.0153 - 1.0)
     38  REAL :: retv
     39  PARAMETER (retv=28.9644/18.0153 - 1.0)
    4040
    41       real zqsat
    42       integer ip
     41  real :: zqsat
     42  integer :: ip
    4343
    44 !    ------------------------------------------------------------------
     44  ! ------------------------------------------------------------------
    4545
    4646
    47       do ip=1,np
     47  do ip=1,np
    4848
    49 !      write(*,*)'kelvin,millibar=',kelvin,millibar
    50 !       write(*,*)'temp,pres=',temp(ip),pres(ip)
     49   ! write(*,*)'kelvin,millibar=',kelvin,millibar
     50   !  write(*,*)'temp,pres=',temp(ip),pres(ip)
    5151
    52          IF (temp(ip) <= rtt) THEN
    53             r3es = r3ies
    54             r4es = r4ies
    55          ELSE
    56             r3es = r3les
    57             r4es = r4les
    58          ENDIF
     52     IF (temp(ip) <= rtt) THEN
     53        r3es = r3ies
     54        r4es = r4ies
     55     ELSE
     56        r3es = r3les
     57        r4es = r4les
     58     ENDIF
    5959
    60          zqsat=r2es/pres(ip)*EXP(r3es*(temp(ip)-rtt)/(temp(ip)-r4es))
    61          zqsat=MIN(0.5,ZQSAT)
    62          zqsat=zqsat/(1.-retv *zqsat)
     60     zqsat=r2es/pres(ip)*EXP(r3es*(temp(ip)-rtt)/(temp(ip)-r4es))
     61     zqsat=MIN(0.5,ZQSAT)
     62     zqsat=zqsat/(1.-retv *zqsat)
    6363
    64          qsat(ip)= zqsat
    65 !      write(*,*)'qsat=',qsat(ip)
     64     qsat(ip)= zqsat
     65   ! write(*,*)'qsat=',qsat(ip)
    6666
    67       enddo
     67  enddo
    6868
    69       RETURN
    70       END
     69
     70END SUBROUTINE q_sat

  • LMDZ6/branches/Amaury_dev/libf/misc/ran1.f90

    r5104 r5105  
    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<0.OR.IFF==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
    24           IX1=MOD(IA1*IX1+IC1,M1)
    25           IX2=MOD(IA2*IX2+IC2,M2)
    26           R(J)=(REAL(IX1)+REAL(IX2)*RM2)*RM1
    27       END DO
    28         IDUM=1
    29       ENDIF
     16  IF (IDUM<0.OR.IFF==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>97.OR.J<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>97.OR.J<1) stop 1
     35  RAN1=R(J)
     36  R(J)=(REAL(IX1)+REAL(IX2)*RM2)*RM1
     37
     38END FUNCTION RAN1

  • LMDZ6/branches/Amaury_dev/libf/misc/sort.f90

    r5104 r5105  
    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)<=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)<=p) THEN
     25       k=j
     26       p=d(j)
     27     ENDIF
     28    ENDDO
    2929
    30        IF(k/=i) THEN
    31          d(k)=d(i)
    32          d(i)=p
    33        ENDIF
    34       ENDDO
     30   IF(k/=i) THEN
     31     d(k)=d(i)
     32     d(i)=p
     33   ENDIF
     34  ENDDO
    3535
    36        RETURN
    37        END
     36
     37END SUBROUTINE sort

  • LMDZ6/branches/Amaury_dev/libf/misc/xercnt.f90

    r5104 r5105  
    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
     62END SUBROUTINE XERCNT

  • LMDZ6/branches/Amaury_dev/libf/misc/xerhlt.f90

    r5104 r5105  
    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/branches/Amaury_dev/libf/misc/xermsg.f90

    r5104 r5105  
    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<-9999999 .OR. NERR>99999999 .OR. NERR==0 .OR.
    208      *   LEVEL<-1 .OR. LEVEL>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==-1 .AND. KOUNT>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<2 .AND. LKNTRL==0) GO TO 30
    242       IF (LEVEL==0 .AND. KOUNT>MAXMES) GO TO 30
    243       IF (LEVEL==1 .AND. KOUNT>MAXMES .AND. MKNTRL==1) GO TO 30
    244       IF (LEVEL==2 .AND. KOUNT>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 /= 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 > 0) THEN
    284 C
    285 C       THE FIRST PART OF THE MESSAGE TELLS ABOUT THE LEVEL.
    286 C
    287          IF (LEVEL <= 0) THEN
    288             TEMP(1:20) = 'INFORMATIVE MESSAGE,'
    289             LTEMP = 20
    290          ELSEIF (LEVEL == 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==2 .AND. LEVEL>=1) .OR.
    301      *       (MKNTRL==1 .AND. LEVEL==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 > 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 > 0) THEN
    329          WRITE (TEMP, '(''ERROR NUMBER = '', I8)') NERR
    330          DO I=16,22
    331             IF (TEMP(I:I) /= ' ') GO TO 20
    332       END DO
    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 /= 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<=0 .OR. (LEVEL==1 .AND. MKNTRL<=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>0 .AND. KOUNT<MAX(1,MAXMES)) THEN
    356          IF (LEVEL == 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<-9999999 .OR. NERR>99999999 .OR. NERR==0 .OR. &
     208        LEVEL<-1 .OR. LEVEL>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==-1 .AND. KOUNT>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<2 .AND. LKNTRL==0) GO TO 30
     242  IF (LEVEL==0 .AND. KOUNT>MAXMES) GO TO 30
     243  IF (LEVEL==1 .AND. KOUNT>MAXMES .AND. MKNTRL==1) GO TO 30
     244  IF (LEVEL==2 .AND. KOUNT>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 /= 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 > 0) THEN
     284  !
     285  !   THE FIRST PART OF THE MESSAGE TELLS ABOUT THE LEVEL.
     286  !
     287     IF (LEVEL <= 0) THEN
     288        TEMP(1:20) = 'INFORMATIVE MESSAGE,'
     289        LTEMP = 20
     290     ELSEIF (LEVEL == 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==2 .AND. LEVEL>=1) .OR. &
     301           (MKNTRL==1 .AND. LEVEL==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 > 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 > 0) THEN
     329     WRITE (TEMP, '(''ERROR NUMBER = '', I8)') NERR
     330     DO I=16,22
     331        IF (TEMP(I:I) /= ' ') 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 /= 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<=0 .OR. (LEVEL==1 .AND. MKNTRL<=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>0 .AND. KOUNT<MAX(1,MAXMES)) THEN
     356     IF (LEVEL == 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
     368END SUBROUTINE XERMSG

  • LMDZ6/branches/Amaury_dev/libf/misc/xerprn.f90

    r5104 r5105  
    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 I=1,NUNIT
    95          IF (IU(I) == 0) IU(I) = N
    96       END DO
    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 < 0 ) THEN
    103          LPREF = LEN(PREFIX)
    104       ELSE
    105          LPREF = NPREF
    106       ENDIF
    107       LPREF = MIN(16, LPREF)
    108       IF (LPREF /= 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 I=1,N
    120          IF (MESSG(LENMSG:LENMSG) /= ' ') GO TO 30
    121          LENMSG = LENMSG - 1
    122       END DO
    123    30 CONTINUE
    124 C
    125 C       IF THE MESSAGE IS ALL BLANKS, THEN PRINT ONE BLANK LINE.
    126 C
    127       IF (LENMSG == 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 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 == 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 < LENMSG+1-NEXTC) THEN
    181             DO I=LPIECE+1,2,-1
    182                IF (MESSG(NEXTC+I-1:NEXTC+I-1) == ' ') 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 == 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 > LWRAP+1) THEN
    199 C
    200 C       LPIECE SHOULD BE SET DOWN TO LWRAP.
    201 C
    202          IDELTA = 0
    203          LPIECE = LWRAP
    204          DO I=LPIECE+1,2,-1
    205             IF (MESSG(NEXTC+I-1:NEXTC+I-1) == ' ') 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 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 I=1,NUNIT
    226          WRITE(IU(I), '(A)') CBUFF(1:LPREF+LPIECE)
    227       END DO
    228 C
    229       IF (NEXTC <= 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) == 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 < 0 ) THEN
     103     LPREF = LEN(PREFIX)
     104  ELSE
     105     LPREF = NPREF
     106  ENDIF
     107  LPREF = MIN(16, LPREF)
     108  IF (LPREF /= 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) /= ' ') 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 == 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 == 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 < LENMSG+1-NEXTC) THEN
     181        DO I=LPIECE+1,2,-1
     182           IF (MESSG(NEXTC+I-1:NEXTC+I-1) == ' ') 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 == 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 > 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) == ' ') 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 <= LENMSG) GO TO 50
     230
     231END SUBROUTINE XERPRN

  • LMDZ6/branches/Amaury_dev/libf/misc/xersve.f90

    r5104 r5105  
    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<=0) THEN
    75 C
    76 C        Dump the table.
    77 C
    78          IF (NMSG==0) RETURN
    79 C
    80 C        Print to each unit.
    81 C
    82          CALL XGETUA (LUN, NUNIT)
    83          DO KUNIT = 1,NUNIT
    84             IUNIT = LUN(KUNIT)
    85             IF (IUNIT==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 I = 1,NMSG
    94                WRITE (IUNIT,9010) LIBTAB(I), SUBTAB(I), MESTAB(I),
    95      *            NERTAB(I),LEVTAB(I),KOUNT(I)
    96       END DO
    97 C
    98 C           Print number of other errors.
    99 C
    100             IF (KOUNTX/=0) WRITE (IUNIT,9020) KOUNTX
    101             WRITE (IUNIT,9030)
    102       END DO
    103 C
    104 C        Clear the error tables.
    105 C
    106          IF (KFLAG==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 I = 1,NMSG
    120             IF (LIB==LIBTAB(I) .AND. SUB==SUBTAB(I) .AND.
    121      *         MES==MESTAB(I) .AND. NERR==NERTAB(I) .AND.
    122      *         LEVEL==LEVTAB(I)) THEN
    123                   KOUNT(I) = KOUNT(I) + 1
    124                   ICOUNT = KOUNT(I)
    125                   RETURN
    126             ENDIF
    127       END DO
    128 C
    129          IF (NMSG<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<=0) THEN
     75  !
     76  !    Dump the table.
     77  !
     78     IF (NMSG==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==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/=0) WRITE (IUNIT,9020) KOUNTX
     101        WRITE (IUNIT,9030)
     102  END DO
     103  !
     104  !    Clear the error tables.
     105  !
     106     IF (KFLAG==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==LIBTAB(I) .AND. SUB==SUBTAB(I) .AND. &
     121              MES==MESTAB(I) .AND. NERR==NERTAB(I) .AND. &
     122              LEVEL==LEVTAB(I)) THEN
     123              KOUNT(I) = KOUNT(I) + 1
     124              ICOUNT = KOUNT(I)
     125              RETURN
     126        ENDIF
     127  END DO
     128  !
     129     IF (NMSG<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/branches/Amaury_dev/libf/misc/xgetua.f90

    r5104 r5105  
    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 I=1,N
    48          INDEX = I+4
    49          IF (I==1) INDEX = 3
    50          IUNITA(I) = J4SAVE(INDEX,0,.FALSE.)
    51       END DO
    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==1) INDEX = 3
     50     IUNITA(I) = J4SAVE(INDEX,0,.FALSE.)
     51  END DO
     52
     53END SUBROUTINE XGETUA
Note: See TracChangeset for help on using the changeset viewer.