Ignore:
Timestamp:
Jul 24, 2024, 1:17:08 PM (12 months ago)
Author:
abarral
Message:

Rename modules in misc from *_mod > lmdz_*
Put cbrt.f90, ch*.f90, pch*.f90 in new lmdz_libmath_pch.f90

Location:
LMDZ6/branches/Amaury_dev/libf/misc
Files:
6 deleted
8 edited
4 moved

Legend:

Unmodified
Added
Removed

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

    r5105 r5113  
    33
    44SUBROUTINE formcoord(unit,n,x,a,rev,text)
    5   implicit none
     5  IMPLICIT NONE
    66  integer :: n,unit,ndec
    77  logical :: rev

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

    r5112 r5113  
    1 MODULE arth_m
     1MODULE lmdz_arth
    22
    3   IMPLICIT NONE
     3  IMPLICIT NONE; PRIVATE
     4  PUBLIC arth
    45
    5   INTEGER, PARAMETER, private:: NPAR_ARTH=16, NPAR2_ARTH=8
     6  INTEGER, PARAMETER :: NPAR_ARTH = 16, NPAR2_ARTH = 8
    67
    78  INTERFACE arth
    8      ! Returns an arithmetic progression, given a first term "first", an
    9      ! increment and a number of terms "n" (including "first").
     9    ! Returns an arithmetic progression, given a first term "first", an
     10    ! increment and a number of terms "n" (including "first").
    1011
    11      MODULE PROCEDURE arth_r, arth_i
    12      ! The difference between the procedures is the kind and type of
    13      ! arguments first and increment and of function result.
     12    MODULE PROCEDURE arth_r, arth_i
     13    ! The difference between the procedures is the kind and type of
     14    ! arguments first and increment and of function result.
    1415  END INTERFACE
    15 
    16   private arth_r, arth_i
    1716
    1817CONTAINS
    1918
    20   pure FUNCTION arth_r(first,increment,n)
     19  PURE FUNCTION arth_r(first, increment, n)
    2120
    22     REAL, INTENT(IN) :: first,increment
     21    REAL, INTENT(IN) :: first, increment
    2322    INTEGER, INTENT(IN) :: n
    2423    REAL arth_r(n)
    2524
    2625    ! Local:
    27     INTEGER :: k,k2
     26    INTEGER :: k, k2
    2827    REAL :: temp
    2928
    3029    !---------------------------------------
    3130
    32     if (n > 0) arth_r(1)=first
    33     if (n <= NPAR_ARTH) then
    34        do k=2,n
    35           arth_r(k)=arth_r(k-1)+increment
    36        END DO
    37     else
    38        do k=2,NPAR2_ARTH
    39           arth_r(k)=arth_r(k-1)+increment
    40        END DO
    41        temp=increment*NPAR2_ARTH
    42        k=NPAR2_ARTH
    43        do
    44           if (k >= n) exit
    45           k2=k+k
    46           arth_r(k+1:min(k2,n)) = temp + arth_r(1:min(k,n-k))
    47           temp=temp+temp
    48           k=k2
    49        END DO
    50     end if
     31    IF (n > 0) arth_r(1) = first
     32    IF (n <= NPAR_ARTH) THEN
     33      DO k = 2, n
     34        arth_r(k) = arth_r(k - 1) + increment
     35      END DO
     36    ELSE
     37      DO k = 2, NPAR2_ARTH
     38        arth_r(k) = arth_r(k - 1) + increment
     39      END DO
     40      temp = increment * NPAR2_ARTH
     41      k = NPAR2_ARTH
     42      DO
     43        IF (k >= n) EXIT
     44        k2 = k + k
     45        arth_r(k + 1:min(k2, n)) = temp + arth_r(1:min(k, n - k))
     46        temp = temp + temp
     47        k = k2
     48      END DO
     49    END IF
    5150
    5251  END FUNCTION arth_r
    … …  
    5453  !*************************************
    5554
    56   pure FUNCTION arth_i(first,increment,n)
     55  PURE FUNCTION arth_i(first, increment, n)
    5756
    58     INTEGER, INTENT(IN) :: first,increment,n
     57    INTEGER, INTENT(IN) :: first, increment, n
    5958    INTEGER arth_i(n)
    6059
    6160    ! Local:
    62     INTEGER :: k,k2,temp
     61    INTEGER :: k, k2, temp
    6362
    6463    !---------------------------------------
    6564
    66     if (n > 0) arth_i(1)=first
    67     if (n <= NPAR_ARTH) then
    68        do k=2,n
    69           arth_i(k)=arth_i(k-1)+increment
    70        END DO
    71     else
    72        do k=2,NPAR2_ARTH
    73           arth_i(k)=arth_i(k-1)+increment
    74        END DO
    75        temp=increment*NPAR2_ARTH
    76        k=NPAR2_ARTH
    77        do
    78           if (k >= n) exit
    79           k2=k+k
    80           arth_i(k+1:min(k2,n))=temp+arth_i(1:min(k,n-k))
    81           temp=temp+temp
    82           k=k2
    83        END DO
    84     end if
     65    IF (n > 0) arth_i(1) = first
     66    IF (n <= NPAR_ARTH) THEN
     67      DO k = 2, n
     68        arth_i(k) = arth_i(k - 1) + increment
     69      END DO
     70    ELSE
     71      DO k = 2, NPAR2_ARTH
     72        arth_i(k) = arth_i(k - 1) + increment
     73      END DO
     74      temp = increment * NPAR2_ARTH
     75      k = NPAR2_ARTH
     76      DO
     77        IF (k >= n) EXIT
     78        k2 = k + k
     79        arth_i(k + 1:min(k2, n)) = temp + arth_i(1:min(k, n - k))
     80        temp = temp + temp
     81        k = k2
     82      END DO
     83    END IF
    8584
    8685  END FUNCTION arth_i
    8786
    88 END MODULE arth_m
     87END MODULE lmdz_arth

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

    r5112 r5113  
    11! $Id$
    2 MODULE assert_m
     2MODULE lmdz_assert
    33
    4   implicit none
     4  IMPLICIT NONE
    55
    66  INTERFACE assert
    7      MODULE PROCEDURE assert1,assert2,assert3,assert4,assert_v
     7    MODULE PROCEDURE assert1, assert2, assert3, assert4, assert_v
    88  END INTERFACE
    99
    10   private assert1,assert2,assert3,assert4,assert_v
     10  PRIVATE assert1, assert2, assert3, assert4, assert_v
    1111
    1212CONTAINS
    1313
    14   SUBROUTINE assert1(n1,string)
    15     CHARACTER(LEN=*), INTENT(IN) :: string
     14  SUBROUTINE assert1(n1, string)
     15    CHARACTER(LEN = *), INTENT(IN) :: string
    1616    LOGICAL, INTENT(IN) :: n1
    1717    if (.not. n1) then
    18        write (*,*) 'nrerror: an assertion failed with this tag:', &
    19             string
    20        print *, 'program terminated by assert1'
    21        stop 1
     18      write (*, *) 'nrerror: an assertion failed with this tag:', &
     19              string
     20      print *, 'program terminated by assert1'
     21      stop 1
    2222    end if
    2323  END SUBROUTINE assert1
    2424  !BL
    25   SUBROUTINE assert2(n1,n2,string)
    26     CHARACTER(LEN=*), INTENT(IN) :: string
    27     LOGICAL, INTENT(IN) :: n1,n2
     25  SUBROUTINE assert2(n1, n2, string)
     26    CHARACTER(LEN = *), INTENT(IN) :: string
     27    LOGICAL, INTENT(IN) :: n1, n2
    2828    if (.not. (n1 .and. n2)) then
    29        write (*,*) 'nrerror: an assertion failed with this tag:', &
    30             string
    31        print *, 'program terminated by assert2'
    32        stop 1
     29      write (*, *) 'nrerror: an assertion failed with this tag:', &
     30              string
     31      print *, 'program terminated by assert2'
     32      stop 1
    3333    end if
    3434  END SUBROUTINE assert2
    3535  !BL
    36   SUBROUTINE assert3(n1,n2,n3,string)
    37     CHARACTER(LEN=*), INTENT(IN) :: string
    38     LOGICAL, INTENT(IN) :: n1,n2,n3
     36  SUBROUTINE assert3(n1, n2, n3, string)
     37    CHARACTER(LEN = *), INTENT(IN) :: string
     38    LOGICAL, INTENT(IN) :: n1, n2, n3
    3939    if (.not. (n1 .and. n2 .and. n3)) then
    40        write (*,*) 'nrerror: an assertion failed with this tag:', &
    41             string
    42        print *, 'program terminated by assert3'
    43        stop 1
     40      write (*, *) 'nrerror: an assertion failed with this tag:', &
     41              string
     42      print *, 'program terminated by assert3'
     43      stop 1
    4444    end if
    4545  END SUBROUTINE assert3
    4646  !BL
    47   SUBROUTINE assert4(n1,n2,n3,n4,string)
    48     CHARACTER(LEN=*), INTENT(IN) :: string
    49     LOGICAL, INTENT(IN) :: n1,n2,n3,n4
     47  SUBROUTINE assert4(n1, n2, n3, n4, string)
     48    CHARACTER(LEN = *), INTENT(IN) :: string
     49    LOGICAL, INTENT(IN) :: n1, n2, n3, n4
    5050    if (.not. (n1 .and. n2 .and. n3 .and. n4)) then
    51        write (*,*) 'nrerror: an assertion failed with this tag:', &
    52             string
    53        print *, 'program terminated by assert4'
    54        stop 1
     51      write (*, *) 'nrerror: an assertion failed with this tag:', &
     52              string
     53      print *, 'program terminated by assert4'
     54      stop 1
    5555    end if
    5656  END SUBROUTINE assert4
    5757  !BL
    58   SUBROUTINE assert_v(n,string)
    59     CHARACTER(LEN=*), INTENT(IN) :: string
     58  SUBROUTINE assert_v(n, string)
     59    CHARACTER(LEN = *), INTENT(IN) :: string
    6060    LOGICAL, DIMENSION(:), INTENT(IN) :: n
    6161    if (.not. all(n)) then
    62        write (*,*) 'nrerror: an assertion failed with this tag:', &
    63             string
    64        print *, 'program terminated by assert_v'
    65        stop 1
     62      write (*, *) 'nrerror: an assertion failed with this tag:', &
     63              string
     64      print *, 'program terminated by assert_v'
     65      stop 1
    6666    end if
    6767  END SUBROUTINE assert_v
    6868
    69 END MODULE assert_m
     69END MODULE lmdz_assert

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

    r5112 r5113  
    1 ! $Id$
    2 MODULE assert_eq_m
     1MODULE lmdz_assert_eq
    32
    4   implicit none
     3  IMPLICIT NONE
    54
    65  INTERFACE assert_eq
    7      MODULE PROCEDURE assert_eq2,assert_eq3,assert_eq4,assert_eqn
     6    MODULE PROCEDURE assert_eq2, assert_eq3, assert_eq4, assert_eqn
    87  END INTERFACE
    98
    10   private assert_eq2,assert_eq3,assert_eq4,assert_eqn
     9  PRIVATE assert_eq2, assert_eq3, assert_eq4, assert_eqn
    1110
    1211CONTAINS
    1312
    14   FUNCTION assert_eq2(n1,n2,string)
    15     CHARACTER(LEN=*), INTENT(IN) :: string
    16     INTEGER, INTENT(IN) :: n1,n2
     13  FUNCTION assert_eq2(n1, n2, string)
     14    CHARACTER(LEN = *), INTENT(IN) :: string
     15    INTEGER, INTENT(IN) :: n1, n2
    1716    INTEGER  assert_eq2
    18     if (n1 == n2) then
    19        assert_eq2=n1
    20     else
    21        write (*,*) 'nrerror: an assert_eq failed with this tag: ', &
    22             string
    23        print *, 'program terminated by assert_eq2'
    24        stop 1
    25     end if
     17    IF (n1 == n2) THEN
     18      assert_eq2 = n1
     19    ELSE
     20      WRITE (*, *) 'nrerror: an assert_eq failed with this tag: ', &
     21              string
     22      PRINT *, 'program terminated by assert_eq2'
     23      STOP 1
     24    END IF
    2625  END FUNCTION assert_eq2
    2726  !BL
    28   FUNCTION assert_eq3(n1,n2,n3,string)
    29     CHARACTER(LEN=*), INTENT(IN) :: string
    30     INTEGER, INTENT(IN) :: n1,n2,n3
     27  FUNCTION assert_eq3(n1, n2, n3, string)
     28    CHARACTER(LEN = *), INTENT(IN) :: string
     29    INTEGER, INTENT(IN) :: n1, n2, n3
    3130    INTEGER  assert_eq3
    32     if (n1 == n2 .and. n2 == n3) then
    33        assert_eq3=n1
    34     else
    35        write (*,*) 'nrerror: an assert_eq failed with this tag: ', &
    36             string
    37        print *, 'program terminated by assert_eq3'
    38        stop 1
    39     end if
     31    IF (n1 == n2 .and. n2 == n3) THEN
     32      assert_eq3 = n1
     33    ELSE
     34      WRITE (*, *) 'nrerror: an assert_eq failed with this tag: ', &
     35              string
     36      PRINT *, 'program terminated by assert_eq3'
     37      STOP 1
     38    END IF
    4039  END FUNCTION assert_eq3
    4140  !BL
    42   FUNCTION assert_eq4(n1,n2,n3,n4,string)
    43     CHARACTER(LEN=*), INTENT(IN) :: string
    44     INTEGER, INTENT(IN) :: n1,n2,n3,n4
     41  FUNCTION assert_eq4(n1, n2, n3, n4, string)
     42    CHARACTER(LEN = *), INTENT(IN) :: string
     43    INTEGER, INTENT(IN) :: n1, n2, n3, n4
    4544    INTEGER  assert_eq4
    46     if (n1 == n2 .and. n2 == n3 .and. n3 == n4) then
    47        assert_eq4=n1
    48     else
    49        write (*,*) 'nrerror: an assert_eq failed with this tag: ', &
    50             string,n1,n2,n3,n4
    51        print *, 'program terminated by assert_eq4'
    52        stop 1
    53     end if
     45    IF (n1 == n2 .and. n2 == n3 .and. n3 == n4) THEN
     46      assert_eq4 = n1
     47    ELSE
     48      WRITE (*, *) 'nrerror: an assert_eq failed with this tag: ', &
     49              string, n1, n2, n3, n4
     50      PRINT *, 'program terminated by assert_eq4'
     51      STOP 1
     52    END IF
    5453  END FUNCTION assert_eq4
    5554  !BL
    56   FUNCTION assert_eqn(nn,string)
    57     CHARACTER(LEN=*), INTENT(IN) :: string
     55  FUNCTION assert_eqn(nn, string)
     56    CHARACTER(LEN = *), INTENT(IN) :: string
    5857    INTEGER, DIMENSION(:), INTENT(IN) :: nn
    5958    INTEGER  assert_eqn
    60     if (all(nn(2:) == nn(1))) then
    61        assert_eqn=nn(1)
    62     else
    63        write (*,*) 'nrerror: an assert_eq failed with this tag: ', &
    64             string
    65        print *, 'program terminated by assert_eqn'
    66        stop 1
    67     end if
     59    IF (all(nn(2:) == nn(1))) THEN
     60      assert_eqn = nn(1)
     61    ELSE
     62      WRITE (*, *) 'nrerror: an assert_eq failed with this tag: ', &
     63              string
     64      PRINT *, 'program terminated by assert_eqn'
     65      STOP 1
     66    END IF
    6867  END FUNCTION assert_eqn
    6968
    70 END MODULE assert_eq_m
     69END MODULE lmdz_assert_eq

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

    r5105 r5113  
    1 FUNCTION cbrt(x)
    2   ! ! Return the cubic root for positive or negative x
    3   IMPLICIT NONE
    4 
    5   REAL :: x,cbrt
    6 
    7   cbrt=sign(1.,x)*(abs(x)**(1./3.))
    8 
    9 
    10 END FUNCTION cbrt
    11 
     1! This module encapsulates misc. math functions
     2
     3MODULE lmdz_libmath_pch
     4  IMPLICIT NONE; PRIVATE
     5  PUBLIC pchfe_95, pchsp_95
     6CONTAINS
     7  REAL FUNCTION cbrt(x)
     8    ! Return the cubic root for positive or negative x
     9    REAL :: x
     10
     11    cbrt = sign(1., x) * (abs(x)**(1. / 3.))
     12  END FUNCTION cbrt
     13
     14  SUBROUTINE CHFEV(X1, X2, F1, F2, D1, D2, NE, XE, FE, NEXT, IERR)
     15    !***BEGIN PROLOGUE  CHFEV
     16    !***PURPOSE  Evaluate a cubic polynomial given in Hermite form at an
     17    ! array of points.  While designed for use by PCHFE, it may
     18    ! be useful directly as an evaluator for a piecewise cubic
     19    ! Hermite function in applications, such as graphing, where
     20    ! the interval is known in advance.
     21    !***LIBRARY   SLATEC (PCHIP)
     22    !***CATEGORY  E3
     23    !***TYPE      SINGLE PRECISION (CHFEV-S, DCHFEV-D)
     24    !***KEYWORDS  CUBIC HERMITE EVALUATION, CUBIC POLYNOMIAL EVALUATION,
     25    !  PCHIP
     26    !***AUTHOR  Fritsch, F. N., (LLNL)
     27    !  Lawrence Livermore National Laboratory
     28    !  P.O. Box 808  (L-316)
     29    !  Livermore, CA  94550
     30    !  FTS 532-4275, (510) 422-4275
     31    !***DESCRIPTION
     32    !
     33    !      CHFEV:  Cubic Hermite Function EValuator
     34    !
     35    ! Evaluates the cubic polynomial determined by function values
     36    ! F1,F2 and derivatives D1,D2 on interval (X1,X2) at the points
     37    ! XE(J), J=1(1)NE.
     38    !
     39    ! ----------------------------------------------------------------------
     40    !
     41    !  Calling sequence:
     42    !
     43    !    INTEGER  NE, NEXT(2), IERR
     44    !    REAL  X1, X2, F1, F2, D1, D2, XE(NE), FE(NE)
     45    !
     46    !    CALL  CHFEV (X1,X2, F1,F2, D1,D2, NE, XE, FE, NEXT, IERR)
     47    !
     48    !   Parameters:
     49    !
     50    ! X1,X2 -- (input) endpoints of interval of definition of cubic.
     51    !       (Error return if  X1.EQ.X2 .)
     52    !
     53    ! F1,F2 -- (input) values of function at X1 and X2, respectively.
     54    !
     55    ! D1,D2 -- (input) values of derivative at X1 and X2, respectively.
     56    !
     57    ! NE -- (input) number of evaluation points.  (Error return if
     58    !       NE.LT.1 .)
     59    !
     60    ! XE -- (input) real array of points at which the function is to be
     61    !       evaluated.  If any of the XE are outside the interval
     62    !       [X1,X2], a warning error is returned in NEXT.
     63    !
     64    ! FE -- (output) real array of values of the cubic function defined
     65    !       by  X1,X2, F1,F2, D1,D2  at the points  XE.
     66    !
     67    ! NEXT -- (output) integer array indicating number of extrapolation
     68    !       points:
     69    !        NEXT(1) = number of evaluation points to left of interval.
     70    !        NEXT(2) = number of evaluation points to right of interval.
     71    !
     72    ! IERR -- (output) error flag.
     73    !       Normal return:
     74    !          IERR = 0  (no errors).
     75    !       "Recoverable" errors:
     76    !          IERR = -1  if NE.LT.1 .
     77    !          IERR = -2  if X1.EQ.X2 .
     78    !            (The FE-array has not been changed in either case.)
     79    !
     80    !***REFERENCES  (NONE)
     81    !***ROUTINES CALLED  XERMSG
     82    !***REVISION HISTORY  (YYMMDD)
     83    !   811019  DATE WRITTEN
     84    !   820803  Minor cosmetic changes for release 1.
     85    !   890411  Added SAVE statements (Vers. 3.2).
     86    !   890531  Changed all specific intrinsics to generic.  (WRB)
     87    !   890703  Corrected category record.  (WRB)
     88    !   890703  REVISION DATE from Version 3.2
     89    !   891214  Prologue converted to Version 4.0 format.  (BAB)
     90    !   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
     91    !***END PROLOGUE  CHFEV
     92    !  Programming notes:
     93    !
     94    ! To produce a double precision version, simply:
     95    !    a. Change CHFEV to DCHFEV wherever it occurs,
     96    !    b. Change the real declaration to double precision, and
     97    !    c. Change the constant ZERO to double precision.
     98    !
     99    !  DECLARE ARGUMENTS.
     100    !
     101    INTEGER :: NE, NEXT(2), IERR
     102    REAL :: X1, X2, F1, F2, D1, D2, XE(*), FE(*)
     103    !
     104    !  DECLARE LOCAL VARIABLES.
     105    !
     106    INTEGER :: I
     107    REAL :: C2, C3, DEL1, DEL2, DELTA, H, X, XMI, XMA, ZERO
     108    SAVE ZERO
     109    DATA  ZERO /0./
     110    !
     111    !  VALIDITY-CHECK ARGUMENTS.
     112    !
     113    !***FIRST EXECUTABLE STATEMENT  CHFEV
     114    IF (NE < 1)  GO TO 5001
     115    H = X2 - X1
     116    IF (H == ZERO)  GO TO 5002
     117    !
     118    !  INITIALIZE.
     119    !
     120    IERR = 0
     121    NEXT(1) = 0
     122    NEXT(2) = 0
     123    XMI = MIN(ZERO, H)
     124    XMA = MAX(ZERO, H)
     125    !
     126    !  COMPUTE CUBIC COEFFICIENTS (EXPANDED ABOUT X1).
     127    !
     128    DELTA = (F2 - F1) / H
     129    DEL1 = (D1 - DELTA) / H
     130    DEL2 = (D2 - DELTA) / H
     131    ! (DELTA IS NO LONGER NEEDED.)
     132    C2 = -(DEL1 + DEL1 + DEL2)
     133    C3 = (DEL1 + DEL2) / H
     134    ! (H, DEL1 AND DEL2 ARE NO LONGER NEEDED.)
     135    !
     136    !  EVALUATION LOOP.
     137    !
     138    DO I = 1, NE
     139      X = XE(I) - X1
     140      FE(I) = F1 + X * (D1 + X * (C2 + X * C3))
     141      ! COUNT EXTRAPOLATION POINTS.
     142      IF (X<XMI)  NEXT(1) = NEXT(1) + 1
     143      IF (X>XMA)  NEXT(2) = NEXT(2) + 1
     144      ! (NOTE REDUNDANCY--IF EITHER CONDITION IS TRUE, OTHER IS FALSE.)
     145    END DO
     146    !
     147    !  NORMAL RETURN.
     148    !
     149    RETURN
     150    !
     151    !  ERROR RETURNS.
     152    !
     153    5001   CONTINUE
     154    ! NE.LT.1 RETURN.
     155    IERR = -1
     156    CALL XERMSG ('SLATEC', 'CHFEV', &
     157            'NUMBER OF EVALUATION POINTS LESS THAN ONE', IERR, 1)
     158    RETURN
     159    !
     160    5002   CONTINUE
     161    ! X1.EQ.X2 RETURN.
     162    IERR = -2
     163    CALL XERMSG ('SLATEC', 'CHFEV', 'INTERVAL ENDPOINTS EQUAL', IERR, &
     164            1)
     165    RETURN
     166    !------------- LAST LINE OF CHFEV FOLLOWS ------------------------------
     167  END SUBROUTINE CHFEV
     168
     169  SUBROUTINE PCHFE(N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)
     170    !***BEGIN PROLOGUE  PCHFE
     171    !***PURPOSE  Evaluate a piecewise cubic Hermite function at an array of
     172    ! points.  May be used by itself for Hermite interpolation,
     173    ! or as an evaluator for PCHIM or PCHIC.
     174    !***LIBRARY   SLATEC (PCHIP)
     175    !***CATEGORY  E3
     176    !***TYPE      SINGLE PRECISION (PCHFE-S, DPCHFE-D)
     177    !***KEYWORDS  CUBIC HERMITE EVALUATION, HERMITE INTERPOLATION, PCHIP,
     178    !  PIECEWISE CUBIC EVALUATION
     179    !***AUTHOR  Fritsch, F. N., (LLNL)
     180    !  Lawrence Livermore National Laboratory
     181    !  P.O. Box 808  (L-316)
     182    !  Livermore, CA  94550
     183    !  FTS 532-4275, (510) 422-4275
     184    !***DESCRIPTION
     185    !
     186    !      PCHFE:  Piecewise Cubic Hermite Function Evaluator
     187    !
     188    ! Evaluates the cubic Hermite function defined by  N, X, F, D  at
     189    ! the points  XE(J), J=1(1)NE.
     190    !
     191    ! To provide compatibility with PCHIM and PCHIC, includes an
     192    ! increment between successive values of the F- and D-arrays.
     193    !
     194    ! ----------------------------------------------------------------------
     195    !
     196    !  Calling sequence:
     197    !
     198    !    PARAMETER  (INCFD = ...)
     199    !    INTEGER  N, NE, IERR
     200    !    REAL  X(N), F(INCFD,N), D(INCFD,N), XE(NE), FE(NE)
     201    !    LOGICAL  SKIP
     202    !
     203    !    CALL  PCHFE (N, X, F, D, INCFD, SKIP, NE, XE, FE, IERR)
     204    !
     205    !   Parameters:
     206    !
     207    ! N -- (input) number of data points.  (Error return if N.LT.2 .)
     208    !
     209    ! X -- (input) real array of independent variable values.  The
     210    !       elements of X must be strictly increasing:
     211    !            X(I-1) .LT. X(I),  I = 2(1)N.
     212    !       (Error return if not.)
     213    !
     214    ! F -- (input) real array of function values.  F(1+(I-1)*INCFD) is
     215    !       the value corresponding to X(I).
     216    !
     217    ! D -- (input) real array of derivative values.  D(1+(I-1)*INCFD) is
     218    !       the value corresponding to X(I).
     219    !
     220    ! INCFD -- (input) increment between successive values in F and D.
     221    !       (Error return if  INCFD.LT.1 .)
     222    !
     223    ! SKIP -- (input/output) logical variable which should be set to
     224    !       .TRUE. if the user wishes to skip checks for validity of
     225    !       preceding parameters, or to .FALSE. otherwise.
     226    !       This will save time in case these checks have already
     227    !       been performed (say, in PCHIM or PCHIC).
     228    !       SKIP will be set to .TRUE. on normal return.
     229    !
     230    ! NE -- (input) number of evaluation points.  (Error return if
     231    !       NE.LT.1 .)
     232    !
     233    ! XE -- (input) real array of points at which the function is to be
     234    !       evaluated.
     235    !
     236    !      NOTES:
     237    !       1. The evaluation will be most efficient if the elements
     238    !          of XE are increasing relative to X;
     239    !          that is,   XE(J) .GE. X(I)
     240    !          implies    XE(K) .GE. X(I),  all K.GE.J .
     241    !       2. If any of the XE are outside the interval [X(1),X(N)],
     242    !          values are extrapolated from the nearest extreme cubic,
     243    !          and a warning error is returned.
     244    !
     245    ! FE -- (output) real array of values of the cubic Hermite function
     246    !       defined by  N, X, F, D  at the points  XE.
     247    !
     248    ! IERR -- (output) error flag.
     249    !       Normal return:
     250    !          IERR = 0  (no errors).
     251    !       Warning error:
     252    !          IERR.GT.0  means that extrapolation was performed at
     253    !             IERR points.
     254    !       "Recoverable" errors:
     255    !          IERR = -1  if N.LT.2 .
     256    !          IERR = -2  if INCFD.LT.1 .
     257    !          IERR = -3  if the X-array is not strictly increasing.
     258    !          IERR = -4  if NE.LT.1 .
     259    !         (The FE-array has not been changed in any of these cases.)
     260    !           NOTE:  The above errors are checked in the order listed,
     261    !               and following arguments have **NOT** been validated.
     262    !
     263    !***REFERENCES  (NONE)
     264    !***ROUTINES CALLED  CHFEV, XERMSG
     265    !***REVISION HISTORY  (YYMMDD)
     266    !   811020  DATE WRITTEN
     267    !   820803  Minor cosmetic changes for release 1.
     268    !   870707  Minor cosmetic changes to prologue.
     269    !   890531  Changed all specific intrinsics to generic.  (WRB)
     270    !   890831  Modified array declarations.  (WRB)
     271    !   890831  REVISION DATE from Version 3.2
     272    !   891214  Prologue converted to Version 4.0 format.  (BAB)
     273    !   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
     274    !***END PROLOGUE  PCHFE
     275    !  Programming notes:
     276    !
     277    ! 1. To produce a double precision version, simply:
     278    !    a. Change PCHFE to DPCHFE, and CHFEV to DCHFEV, wherever they
     279    !       occur,
     280    !    b. Change the real declaration to double precision,
     281    !
     282    ! 2. Most of the coding between the CALL to CHFEV and the end of
     283    !    the IR-loop could be eliminated if it were permissible to
     284    !    assume that XE is ordered relative to X.
     285    !
     286    ! 3. CHFEV does not assume that X1 is less than X2.  thus, it would
     287    !    be possible to write a version of PCHFE that assumes a strict-
     288    !    ly decreasing X-array by simply running the IR-loop backwards
     289    !    (and reversing the order of appropriate tests).
     290    !
     291    ! 4. The present code has a minor bug, which I have decided is not
     292    !    worth the effort that would be required to fix it.
     293    !    If XE contains points in [X(N-1),X(N)], followed by points .LT.
     294    !    X(N-1), followed by points .GT.X(N), the extrapolation points
     295    !    will be counted (at least) twice in the total returned in IERR.
     296    !
     297    !  DECLARE ARGUMENTS.
     298    !
     299    INTEGER :: N, INCFD, NE, IERR
     300    REAL :: X(*), F(INCFD, *), D(INCFD, *), XE(*), FE(*)
     301    LOGICAL :: SKIP
     302    !
     303    !  DECLARE LOCAL VARIABLES.
     304    !
     305    INTEGER :: I, IERC, IR, J, JFIRST, NEXT(2), NJ
     306    !
     307    !  VALIDITY-CHECK ARGUMENTS.
     308    !
     309    !***FIRST EXECUTABLE STATEMENT  PCHFE
     310    IF (SKIP)  GO TO 5
     311    !
     312    IF (N<2)  GO TO 5001
     313    IF (INCFD<1)  GO TO 5002
     314    DO I = 2, N
     315      IF (X(I)<=X(I - 1))  GO TO 5003
     316    END DO
     317    !
     318    !  FUNCTION DEFINITION IS OK, GO ON.
     319    !
     320    5   CONTINUE
     321    IF (NE<1)  GO TO 5004
     322    IERR = 0
     323    SKIP = .TRUE.
     324    !
     325    !  LOOP OVER INTERVALS.        (   INTERVAL INDEX IS  IL = IR-1  . )
     326    !                              ( INTERVAL IS X(IL).LE.X.LT.X(IR) . )
     327    JFIRST = 1
     328    IR = 2
     329    10   CONTINUE
     330    !
     331    ! SKIP OUT OF LOOP IF HAVE PROCESSED ALL EVALUATION POINTS.
     332    !
     333    IF (JFIRST > NE)  GO TO 5000
     334    !
     335    ! LOCATE ALL POINTS IN INTERVAL.
     336    !
     337    DO J = JFIRST, NE
     338      IF (XE(J) >= X(IR))  GO TO 30
     339    END DO
     340    J = NE + 1
     341    GO TO 40
     342    !
     343    ! HAVE LOCATED FIRST POINT BEYOND INTERVAL.
     344    !
     345    30   CONTINUE
     346    IF (IR == N)  J = NE + 1
     347    !
     348    40   CONTINUE
     349    NJ = J - JFIRST
     350    !
     351    ! SKIP EVALUATION IF NO POINTS IN INTERVAL.
     352    !
     353    IF (NJ == 0)  GO TO 50
     354    !
     355    ! EVALUATE CUBIC AT XE(I),  I = JFIRST (1) J-1 .
     356    !
     357    !   ----------------------------------------------------------------
     358    CALL CHFEV (X(IR - 1), X(IR), F(1, IR - 1), F(1, IR), D(1, IR - 1), D(1, IR), &
     359            NJ, XE(JFIRST), FE(JFIRST), NEXT, IERC)
     360    ! ----------------------------------------------------------------
     361    IF (IERC < 0)  GO TO 5005
     362    !
     363    IF (NEXT(2) == 0)  GO TO 42
     364    ! IF (NEXT(2) .GT. 0)  THEN
     365    !    IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(2) TO THE
     366    !    RIGHT OF X(IR).
     367    !
     368    IF (IR < N)  GO TO 41
     369    ! IF (IR .EQ. N)  THEN
     370    !    THESE ARE ACTUALLY EXTRAPOLATION POINTS.
     371    IERR = IERR + NEXT(2)
     372    GO TO 42
     373    41    CONTINUE
     374    ! ELSE
     375    !    WE SHOULD NEVER HAVE GOTTEN HERE.
     376    GO TO 5005
     377    ! ENDIF
     378    ! ENDIF
     379    42   CONTINUE
     380    !
     381    IF (NEXT(1) == 0)  GO TO 49
     382    ! IF (NEXT(1) .GT. 0)  THEN
     383    !    IN THE CURRENT SET OF XE-POINTS, THERE ARE NEXT(1) TO THE
     384    !    LEFT OF X(IR-1).
     385    !
     386    IF (IR > 2)  GO TO 43
     387    ! IF (IR .EQ. 2)  THEN
     388    !    THESE ARE ACTUALLY EXTRAPOLATION POINTS.
     389    IERR = IERR + NEXT(1)
     390    GO TO 49
     391    43    CONTINUE
     392    ! ELSE
     393    !    XE IS NOT ORDERED RELATIVE TO X, SO MUST ADJUST
     394    !    EVALUATION INTERVAL.
     395    !
     396    !          FIRST, LOCATE FIRST POINT TO LEFT OF X(IR-1).
     397    DO I = JFIRST, J - 1
     398      IF (XE(I) < X(IR - 1))  GO TO 45
     399    END DO
     400    ! NOTE-- CANNOT DROP THROUGH HERE UNLESS THERE IS AN ERROR
     401    !        IN CHFEV.
     402    GO TO 5005
     403    !
     404    45       CONTINUE
     405    ! RESET J.  (THIS WILL BE THE NEW JFIRST.)
     406    J = I
     407    !
     408    !          NOW FIND OUT HOW FAR TO BACK UP IN THE X-ARRAY.
     409    DO I = 1, IR - 1
     410      IF (XE(J) < X(I)) GO TO 47
     411    END DO
     412    ! NB-- CAN NEVER DROP THROUGH HERE, SINCE XE(J).LT.X(IR-1).
     413    !
     414    47       CONTINUE
     415    ! AT THIS POINT, EITHER  XE(J) .LT. X(1)
     416    !    OR      X(I-1) .LE. XE(J) .LT. X(I) .
     417    ! RESET IR, RECOGNIZING THAT IT WILL BE INCREMENTED BEFORE
     418    ! CYCLING.
     419    IR = MAX(1, I - 1)
     420    ! ENDIF
     421    ! ENDIF
     422    49   CONTINUE
     423    !
     424    JFIRST = J
     425    !
     426    ! END OF IR-LOOP.
     427    !
     428    50   CONTINUE
     429    IR = IR + 1
     430    IF (IR <= N)  GO TO 10
     431    !
     432    !  NORMAL RETURN.
     433    !
     434    5000   CONTINUE
     435    RETURN
     436    !
     437    !  ERROR RETURNS.
     438    !
     439    5001   CONTINUE
     440    ! N.LT.2 RETURN.
     441    IERR = -1
     442    CALL XERMSG ('SLATEC', 'PCHFE', &
     443            'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
     444    RETURN
     445    !
     446    5002   CONTINUE
     447    ! INCFD.LT.1 RETURN.
     448    IERR = -2
     449    CALL XERMSG ('SLATEC', 'PCHFE', 'INCREMENT LESS THAN ONE', IERR, &
     450            1)
     451    RETURN
     452    !
     453    5003   CONTINUE
     454    ! X-ARRAY NOT STRICTLY INCREASING.
     455    IERR = -3
     456    CALL XERMSG ('SLATEC', 'PCHFE', 'X-ARRAY NOT STRICTLY INCREASING' &
     457            , IERR, 1)
     458    RETURN
     459    !
     460    5004   CONTINUE
     461    ! NE.LT.1 RETURN.
     462    IERR = -4
     463    CALL XERMSG ('SLATEC', 'PCHFE', &
     464            'NUMBER OF EVALUATION POINTS LESS THAN ONE', IERR, 1)
     465    RETURN
     466    !
     467    5005   CONTINUE
     468    ! ERROR RETURN FROM CHFEV.
     469    !   *** THIS CASE SHOULD NEVER OCCUR ***
     470    IERR = -5
     471    CALL XERMSG ('SLATEC', 'PCHFE', &
     472            'ERROR RETURN FROM CHFEV -- FATAL', IERR, 2)
     473    RETURN
     474    !------------- LAST LINE OF PCHFE FOLLOWS ------------------------------
     475  END SUBROUTINE PCHFE
     476
     477  SUBROUTINE PCHFE_95(X, F, D, SKIP, XE, FE, IERR)
     478
     479    ! PURPOSE  Evaluate a piecewise cubic Hermite function at an array of
     480    !            points.  May be used by itself for Hermite interpolation,
     481    !            or as an evaluator for PCHIM or PCHIC.
     482    ! CATEGORY  E3
     483    ! KEYWORDS  CUBIC HERMITE EVALUATION, HERMITE INTERPOLATION, PCHIP,
     484    !             PIECEWISE CUBIC EVALUATION
     485
     486    !          PCHFE:  Piecewise Cubic Hermite Function Evaluator
     487    ! Evaluates the cubic Hermite function defined by  X, F, D  at
     488    ! the points  XE.
     489
     490    use lmdz_assert_eq, only: assert_eq
     491
     492    REAL, intent(in) :: X(:) ! real array of independent variable values
     493    ! The elements of X must be strictly increasing.
     494
     495    REAL, intent(in) :: F(:) ! real array of function values
     496    ! F(I) is the value corresponding to X(I).
     497
     498    REAL, intent(in) :: D(:) ! real array of derivative values
     499    ! D(I) is the value corresponding to X(I).
     500
     501    LOGICAL, intent(inout) :: SKIP
     502    ! request to skip checks for validity of "x"
     503    ! If "skip" is false then "pchfe" will check that size(x) >= 2 and
     504    ! "x" is in strictly ascending order.
     505    ! Setting "skip" to true will save time in case these checks have
     506    ! already been performed (say, in "PCHIM" or "PCHIC").
     507    ! "SKIP" will be set to TRUE on normal return.
     508
     509    real, intent(in) :: XE(:) ! points at which the function is to be evaluated
     510    ! NOTES:
     511    ! 1. The evaluation will be most efficient if the elements of XE
     512    ! are increasing relative to X.
     513    ! That is,   XE(J) .GE. X(I)
     514    ! implies    XE(K) .GE. X(I),  all K.GE.J
     515    ! 2. If any of the XE are outside the interval [X(1),X(N)], values
     516    ! are extrapolated from the nearest extreme cubic, and a warning
     517    ! error is returned.
     518
     519    real, intent(out) :: FE(:) ! values of the cubic Hermite function
     520    ! defined by X, F, D at the points XE
     521
     522    integer, intent(out) :: IERR ! error flag
     523    ! Normal return:
     524    ! IERR = 0  no error
     525    ! Warning error:
     526    ! IERR > 0  means that extrapolation was performed at IERR points
     527    ! "Recoverable" errors:
     528    !              IERR = -1  if N < 2
     529    !              IERR = -3  if the X-array is not strictly increasing
     530    !              IERR = -4  if NE < 1
     531    ! NOTE: The above errors are checked in the order listed, and
     532    ! following arguments have **NOT** been validated.
     533
     534    ! Variables local to the procedure:
     535
     536    INTEGER  N, NE
     537
     538    !---------------------------------------
     539
     540    n = assert_eq(size(x), size(f), size(d), "PCHFE_95 n")
     541    ne = assert_eq(size(xe), size(fe), "PCHFE_95 ne")
     542    CALL PCHFE(N, X, F, D, 1, SKIP, NE, XE, FE, IERR)
     543
     544  END SUBROUTINE  PCHFE_95
     545
     546  FUNCTION pchsp_95(x, f, ibeg, iend, vc_beg, vc_end)
     547
     548    ! PURPOSE: Set derivatives needed to determine the Hermite
     549    ! representation of the cubic spline interpolant to given data,
     550    ! with specified boundary conditions.
     551
     552    ! Part of the "pchip" package.
     553
     554    ! CATEGORY: E1A
     555
     556    ! KEYWORDS: cubic hermite interpolation, piecewise cubic
     557    ! interpolation, spline interpolation
     558
     559    ! DESCRIPTION: "pchsp" stands for "Piecewise Cubic Hermite Spline"
     560    ! Computes the Hermite representation of the cubic spline
     561    ! interpolant to the data given in X and F satisfying the boundary
     562    ! conditions specified by Ibeg, iend, vc_beg and VC_end.
     563
     564    ! The resulting piecewise cubic Hermite function may be evaluated
     565    ! by "pchfe" or "pchfd".
     566
     567    ! NOTE: This is a modified version of C. de Boor's cubic spline
     568    ! routine "cubspl".
     569
     570    ! REFERENCE: Carl de Boor, A Practical Guide to Splines, Springer,
     571    ! 2001, pages 43-47
     572
     573    use lmdz_assert_eq, only: assert_eq
     574
     575    real, intent(in) :: x(:)
     576    ! independent variable values
     577    ! The elements of X must be strictly increasing:
     578    !                X(I-1) < X(I),  I = 2...N.
     579    !           (Error return if not.)
     580    ! (error if size(x) < 2)
     581
     582    real, intent(in) :: f(:)
     583    !     dependent variable values to be interpolated
     584    !  F(I) is value corresponding to X(I).
     585
     586    INTEGER, intent(in) :: ibeg
     587    !     desired boundary condition at beginning of data
     588
     589    !        IBEG = 0  to set pchsp_95(1) so that the third derivative is con-
     590    !              tinuous at X(2).  This is the "not a knot" condition
     591    !              provided by de Boor's cubic spline routine CUBSPL.
     592    !              This is the default boundary condition.
     593    !        IBEG = 1  if first derivative at X(1) is given in VC_BEG.
     594    !        IBEG = 2  if second derivative at X(1) is given in VC_BEG.
     595    !        IBEG = 3  to use the 3-point difference formula for pchsp_95(1).
     596    !              (Reverts to the default boundary condition if size(x) < 3 .)
     597    !        IBEG = 4  to use the 4-point difference formula for pchsp_95(1).
     598    !              (Reverts to the default boundary condition if size(x) < 4 .)
     599
     600    !          NOTES:
     601    !           1. An error return is taken if IBEG is out of range.
     602    !           2. For the "natural" boundary condition, use IBEG=2 and
     603    !              VC_BEG=0.
     604
     605    INTEGER, intent(in) :: iend
     606    !           IC(2) = IEND, desired condition at end of data.
     607    !  IEND may take on the same values as IBEG, but applied to
     608    !  derivative at X(N). In case IEND = 1 or 2, The value is given in VC_END.
     609
     610    !          NOTES:
     611    !           1. An error return is taken if IEND is out of range.
     612    !           2. For the "natural" boundary condition, use IEND=2 and
     613    !              VC_END=0.
     614
     615    REAL, intent(in), optional :: vc_beg
     616    ! desired boundary value, as indicated above.
     617    !           VC_BEG need be set only if IBEG = 1 or 2 .
     618
     619    REAL, intent(in), optional :: vc_end
     620    ! desired boundary value, as indicated above.
     621    !           VC_END need be set only if Iend = 1 or 2 .
     622
     623    real pchsp_95(size(x))
     624    ! derivative values at the data points
     625    !           These values will determine the cubic spline interpolant
     626    !           with the requested boundary conditions.
     627    !           The value corresponding to X(I) is stored in
     628    !                PCHSP_95(I),  I=1...N.
     629
     630    ! LOCAL VARIABLES:
     631    real wk(2, size(x)) ! real array of working storage
     632    INTEGER n ! number of data points
     633    integer ierr, ic(2)
     634    real vc(2)
     635
     636    !-------------------------------------------------------------------
     637
     638    n = assert_eq(size(x), size(f), "pchsp_95 n")
     639    if ((ibeg == 1 .or. ibeg == 2) .and. .not. present(vc_beg)) then
     640      print *, "vc_beg required for IBEG = 1 or 2"
     641      stop 1
     642    end if
     643    if ((iend == 1 .or. iend == 2) .and. .not. present(vc_end)) then
     644      print *, "vc_end required for IEND = 1 or 2"
     645      stop 1
     646    end if
     647    ic = (/ibeg, iend/)
     648    if (present(vc_beg)) vc(1) = vc_beg
     649    if (present(vc_end)) vc(2) = vc_end
     650    CALL PCHSP(IC, VC, N, X, F, pchsp_95, 1, WK, size(WK), IERR)
     651    if (ierr /= 0) stop 1
     652
     653  END FUNCTION pchsp_95
     654
     655  REAL FUNCTION PCHDF(K, X, S, IERR)
     656    !***BEGIN PROLOGUE  PCHDF
     657    !***SUBSIDIARY
     658    !***PURPOSE  Computes divided differences for PCHCE and PCHSP
     659    !***LIBRARY   SLATEC (PCHIP)
     660    !***TYPE      SINGLE PRECISION (PCHDF-S, DPCHDF-D)
     661    !***AUTHOR  Fritsch, F. N., (LLNL)
     662    !***DESCRIPTION
     663    !
     664    !      PCHDF:   PCHIP Finite Difference Formula
     665    !
     666    ! Uses a divided difference formulation to compute a K-point approx-
     667    ! imation to the derivative at X(K) based on the data in X and S.
     668    !
     669    ! Called by  PCHCE  and  PCHSP  to compute 3- and 4-point boundary
     670    ! derivative approximations.
     671    !
     672    ! ----------------------------------------------------------------------
     673    !
     674    ! On input:
     675    !    K      is the order of the desired derivative approximation.
     676    !           K must be at least 3 (error return if not).
     677    !    X      contains the K values of the independent variable.
     678    !           X need not be ordered, but the values **MUST** be
     679    !           distinct.  (Not checked here.)
     680    !    S      contains the associated slope values:
     681    !              S(I) = (F(I+1)-F(I))/(X(I+1)-X(I)), I=1(1)K-1.
     682    !           (Note that S need only be of length K-1.)
     683    !
     684    ! On return:
     685    !    S      will be destroyed.
     686    !    IERR   will be set to -1 if K.LT.2 .
     687    !    PCHDF  will be set to the desired derivative approximation if
     688    !           IERR=0 or to zero if IERR=-1.
     689    !
     690    ! ----------------------------------------------------------------------
     691    !
     692    !***SEE ALSO  PCHCE, PCHSP
     693    !***REFERENCES  Carl de Boor, A Practical Guide to Splines, Springer-
     694    !             Verlag, New York, 1978, pp. 10-16.
     695    !***ROUTINES CALLED  XERMSG
     696    !***REVISION HISTORY  (YYMMDD)
     697    !   820503  DATE WRITTEN
     698    !   820805  Converted to SLATEC library version.
     699    !   870813  Minor cosmetic changes.
     700    !   890411  Added SAVE statements (Vers. 3.2).
     701    !   890411  REVISION DATE from Version 3.2
     702    !   891214  Prologue converted to Version 4.0 format.  (BAB)
     703    !   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
     704    !   900328  Added TYPE section.  (WRB)
     705    !   910408  Updated AUTHOR and DATE WRITTEN sections in prologue.  (WRB)
     706    !   920429  Revised format and order of references.  (WRB,FNF)
     707    !   930503  Improved purpose.  (FNF)
     708    !***END PROLOGUE  PCHDF
     709    !
     710    !**End
     711    !
     712    !  DECLARE ARGUMENTS.
     713    !
     714    INTEGER :: K, IERR
     715    REAL :: X(K), S(K)
     716    !
     717    !  DECLARE LOCAL VARIABLES.
     718    !
     719    INTEGER :: I, J
     720    REAL :: VALUE, ZERO
     721    SAVE ZERO
     722    DATA  ZERO /0./
     723    !
     724    !  CHECK FOR LEGAL VALUE OF K.
     725    !
     726    !***FIRST EXECUTABLE STATEMENT  PCHDF
     727    IF (K < 3)  GO TO 5001
     728    !
     729    !  COMPUTE COEFFICIENTS OF INTERPOLATING POLYNOMIAL.
     730    !
     731    DO J = 2, K - 1
     732      DO I = 1, K - J
     733        S(I) = (S(I + 1) - S(I)) / (X(I + J) - X(I))
     734      END DO
     735    END DO
     736    !
     737    !  EVALUATE DERIVATIVE AT X(K).
     738    !
     739    VALUE = S(1)
     740    DO I = 2, K - 1
     741      VALUE = S(I) + VALUE * (X(K) - X(I))
     742    END DO
     743    !
     744    !  NORMAL RETURN.
     745    !
     746    IERR = 0
     747    PCHDF = VALUE
     748    RETURN
     749    !
     750    !  ERROR RETURN.
     751    !
     752    5001   CONTINUE
     753    ! K.LT.3 RETURN.
     754    IERR = -1
     755    CALL XERMSG ('SLATEC', 'PCHDF', 'K LESS THAN THREE', IERR, 1)
     756    PCHDF = ZERO
     757    RETURN
     758    !------------- LAST LINE OF PCHDF FOLLOWS ------------------------------
     759  END FUNCTION PCHDF
     760
     761  SUBROUTINE PCHSP(IC, VC, N, X, F, D, INCFD, WK, NWK, IERR)
     762    !***BEGIN PROLOGUE  PCHSP
     763    !***PURPOSE  Set derivatives needed to determine the Hermite represen-
     764    ! tation of the cubic spline interpolant to given data, with
     765    ! specified boundary conditions.
     766    !***LIBRARY   SLATEC (PCHIP)
     767    !***CATEGORY  E1A
     768    !***TYPE      SINGLE PRECISION (PCHSP-S, DPCHSP-D)
     769    !***KEYWORDS  CUBIC HERMITE INTERPOLATION, PCHIP,
     770    !  PIECEWISE CUBIC INTERPOLATION, SPLINE INTERPOLATION
     771    !***AUTHOR  Fritsch, F. N., (LLNL)
     772    !  Lawrence Livermore National Laboratory
     773    !  P.O. Box 808  (L-316)
     774    !  Livermore, CA  94550
     775    !  FTS 532-4275, (510) 422-4275
     776    !***DESCRIPTION
     777    !
     778    !      PCHSP:   Piecewise Cubic Hermite Spline
     779    !
     780    ! Computes the Hermite representation of the cubic spline inter-
     781    ! polant to the data given in X and F satisfying the boundary
     782    ! conditions specified by IC and VC.
     783    !
     784    ! To facilitate two-dimensional applications, includes an increment
     785    ! between successive values of the F- and D-arrays.
     786    !
     787    ! The resulting piecewise cubic Hermite function may be evaluated
     788    ! by PCHFE or PCHFD.
     789    !
     790    ! NOTE:  This is a modified version of C. de Boor's cubic spline
     791    !        routine CUBSPL.
     792    !
     793    ! ----------------------------------------------------------------------
     794    !
     795    !  Calling sequence:
     796    !
     797    !    PARAMETER  (INCFD = ...)
     798    !    INTEGER  IC(2), N, NWK, IERR
     799    !    REAL  VC(2), X(N), F(INCFD,N), D(INCFD,N), WK(NWK)
     800    !
     801    !    CALL  PCHSP (IC, VC, N, X, F, D, INCFD, WK, NWK, IERR)
     802    !
     803    !   Parameters:
     804    !
     805    ! IC -- (input) integer array of length 2 specifying desired
     806    !       boundary conditions:
     807    !       IC(1) = IBEG, desired condition at beginning of data.
     808    !       IC(2) = IEND, desired condition at end of data.
     809    !
     810    !       IBEG = 0  to set D(1) so that the third derivative is con-
     811    !          tinuous at X(2).  This is the "not a knot" condition
     812    !          provided by de Boor's cubic spline routine CUBSPL.
     813    !          < This is the default boundary condition. >
     814    !       IBEG = 1  if first derivative at X(1) is given in VC(1).
     815    !       IBEG = 2  if second derivative at X(1) is given in VC(1).
     816    !       IBEG = 3  to use the 3-point difference formula for D(1).
     817    !                 (Reverts to the default b.c. if N.LT.3 .)
     818    !       IBEG = 4  to use the 4-point difference formula for D(1).
     819    !                 (Reverts to the default b.c. if N.LT.4 .)
     820    !      NOTES:
     821    !       1. An error return is taken if IBEG is out of range.
     822    !       2. For the "natural" boundary condition, use IBEG=2 and
     823    !          VC(1)=0.
     824    !
     825    !       IEND may take on the same values as IBEG, but applied to
     826    !       derivative at X(N).  In case IEND = 1 or 2, the value is
     827    !       given in VC(2).
     828    !
     829    !      NOTES:
     830    !       1. An error return is taken if IEND is out of range.
     831    !       2. For the "natural" boundary condition, use IEND=2 and
     832    !          VC(2)=0.
     833    !
     834    ! VC -- (input) real array of length 2 specifying desired boundary
     835    !       values, as indicated above.
     836    !       VC(1) need be set only if IC(1) = 1 or 2 .
     837    !       VC(2) need be set only if IC(2) = 1 or 2 .
     838    !
     839    ! N -- (input) number of data points.  (Error return if N.LT.2 .)
     840    !
     841    ! X -- (input) real array of independent variable values.  The
     842    !       elements of X must be strictly increasing:
     843    !            X(I-1) .LT. X(I),  I = 2(1)N.
     844    !       (Error return if not.)
     845    !
     846    ! F -- (input) real array of dependent variable values to be inter-
     847    !       polated.  F(1+(I-1)*INCFD) is value corresponding to X(I).
     848    !
     849    ! D -- (output) real array of derivative values at the data points.
     850    !       These values will determine the cubic spline interpolant
     851    !       with the requested boundary conditions.
     852    !       The value corresponding to X(I) is stored in
     853    !            D(1+(I-1)*INCFD),  I=1(1)N.
     854    !       No other entries in D are changed.
     855    !
     856    ! INCFD -- (input) increment between successive values in F and D.
     857    !       This argument is provided primarily for 2-D applications.
     858    !       (Error return if  INCFD.LT.1 .)
     859    !
     860    ! WK -- (scratch) real array of working storage.
     861    !
     862    ! NWK -- (input) length of work array.
     863    !       (Error return if NWK.LT.2*N .)
     864    !
     865    ! IERR -- (output) error flag.
     866    !       Normal return:
     867    !          IERR = 0  (no errors).
     868    !       "Recoverable" errors:
     869    !          IERR = -1  if N.LT.2 .
     870    !          IERR = -2  if INCFD.LT.1 .
     871    !          IERR = -3  if the X-array is not strictly increasing.
     872    !          IERR = -4  if IBEG.LT.0 or IBEG.GT.4 .
     873    !          IERR = -5  if IEND.LT.0 of IEND.GT.4 .
     874    !          IERR = -6  if both of the above are true.
     875    !          IERR = -7  if NWK is too small.
     876    !           NOTE:  The above errors are checked in the order listed,
     877    !               and following arguments have **NOT** been validated.
     878    !         (The D-array has not been changed in any of these cases.)
     879    !          IERR = -8  in case of trouble solving the linear system
     880    !                     for the interior derivative values.
     881    !         (The D-array may have been changed in this case.)
     882    !         (             Do **NOT** use it!                )
     883    !
     884    !***REFERENCES  Carl de Boor, A Practical Guide to Splines, Springer-
     885    !             Verlag, New York, 1978, pp. 53-59.
     886    !***ROUTINES CALLED  PCHDF, XERMSG
     887    !***REVISION HISTORY  (YYMMDD)
     888    !   820503  DATE WRITTEN
     889    !   820804  Converted to SLATEC library version.
     890    !   870707  Minor cosmetic changes to prologue.
     891    !   890411  Added SAVE statements (Vers. 3.2).
     892    !   890703  Corrected category record.  (WRB)
     893    !   890831  Modified array declarations.  (WRB)
     894    !   890831  REVISION DATE from Version 3.2
     895    !   891214  Prologue converted to Version 4.0 format.  (BAB)
     896    !   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
     897    !   920429  Revised format and order of references.  (WRB,FNF)
     898    !***END PROLOGUE  PCHSP
     899    !  Programming notes:
     900    !
     901    ! To produce a double precision version, simply:
     902    !    a. Change PCHSP to DPCHSP wherever it occurs,
     903    !    b. Change the real declarations to double precision, and
     904    !    c. Change the constants ZERO, HALF, ... to double precision.
     905    !
     906    !  DECLARE ARGUMENTS.
     907    !
     908    INTEGER :: IC(2), N, INCFD, NWK, IERR
     909    REAL :: VC(2), X(*), F(INCFD, *), D(INCFD, *), WK(2, *)
     910    !
     911    !  DECLARE LOCAL VARIABLES.
     912    !
     913    INTEGER :: IBEG, IEND, INDEX, J, NM1
     914    REAL :: G, HALF, ONE, STEMP(3), THREE, TWO, XTEMP(4), ZERO
     915    SAVE ZERO, HALF, ONE, TWO, THREE
     916    REAL :: PCHDF
     917    !
     918    DATA  ZERO /0./, HALF /0.5/, ONE /1./, TWO /2./, THREE /3./
     919    !
     920    !  VALIDITY-CHECK ARGUMENTS.
     921    !
     922    !***FIRST EXECUTABLE STATEMENT  PCHSP
     923    IF (N<2)  GO TO 5001
     924    IF (INCFD<1)  GO TO 5002
     925    DO J = 2, N
     926      IF (X(J)<=X(J - 1))  GO TO 5003
     927    END DO
     928    !
     929    IBEG = IC(1)
     930    IEND = IC(2)
     931    IERR = 0
     932    IF ((IBEG<0).OR.(IBEG>4))  IERR = IERR - 1
     933    IF ((IEND<0).OR.(IEND>4))  IERR = IERR - 2
     934    IF (IERR<0)  GO TO 5004
     935    !
     936    !  FUNCTION DEFINITION IS OK -- GO ON.
     937    !
     938    IF (NWK < 2 * N)  GO TO 5007
     939    !
     940    !  COMPUTE FIRST DIFFERENCES OF X SEQUENCE AND STORE IN WK(1,.). ALSO,
     941    !  COMPUTE FIRST DIVIDED DIFFERENCE OF DATA AND STORE IN WK(2,.).
     942    DO J = 2, N
     943      WK(1, J) = X(J) - X(J - 1)
     944      WK(2, J) = (F(1, J) - F(1, J - 1)) / WK(1, J)
     945    END DO
     946    !
     947    !  SET TO DEFAULT BOUNDARY CONDITIONS IF N IS TOO SMALL.
     948    !
     949    IF (IBEG>N)  IBEG = 0
     950    IF (IEND>N)  IEND = 0
     951    !
     952    !  SET UP FOR BOUNDARY CONDITIONS.
     953    !
     954    IF ((IBEG==1).OR.(IBEG==2))  THEN
     955      D(1, 1) = VC(1)
     956    ELSE IF (IBEG > 2)  THEN
     957      ! PICK UP FIRST IBEG POINTS, IN REVERSE ORDER.
     958      DO J = 1, IBEG
     959        INDEX = IBEG - J + 1
     960        ! INDEX RUNS FROM IBEG DOWN TO 1.
     961        XTEMP(J) = X(INDEX)
     962        IF (J < IBEG)  STEMP(J) = WK(2, INDEX)
     963      END DO
     964      ! --------------------------------
     965      D(1, 1) = PCHDF (IBEG, XTEMP, STEMP, IERR)
     966      ! --------------------------------
     967      IF (IERR /= 0)  GO TO 5009
     968      IBEG = 1
     969    ENDIF
     970    !
     971    IF ((IEND==1).OR.(IEND==2))  THEN
     972      D(1, N) = VC(2)
     973    ELSE IF (IEND > 2)  THEN
     974      ! PICK UP LAST IEND POINTS.
     975      DO J = 1, IEND
     976        INDEX = N - IEND + J
     977        ! INDEX RUNS FROM N+1-IEND UP TO N.
     978        XTEMP(J) = X(INDEX)
     979        IF (J < IEND)  STEMP(J) = WK(2, INDEX + 1)
     980      END DO
     981      ! --------------------------------
     982      D(1, N) = PCHDF (IEND, XTEMP, STEMP, IERR)
     983      ! --------------------------------
     984      IF (IERR /= 0)  GO TO 5009
     985      IEND = 1
     986    ENDIF
     987    !
     988    ! --------------------( BEGIN CODING FROM CUBSPL )--------------------
     989    !
     990    !  **** A TRIDIAGONAL LINEAR SYSTEM FOR THE UNKNOWN SLOPES S(J) OF
     991    !  F  AT X(J), J=1,...,N, IS GENERATED AND THEN SOLVED BY GAUSS ELIM-
     992    !  INATION, WITH S(J) ENDING UP IN D(1,J), ALL J.
     993    ! WK(1,.) AND WK(2,.) ARE USED FOR TEMPORARY STORAGE.
     994    !
     995    !  CONSTRUCT FIRST EQUATION FROM FIRST BOUNDARY CONDITION, OF THE FORM
     996    !         WK(2,1)*S(1) + WK(1,1)*S(2) = D(1,1)
     997    !
     998    IF (IBEG == 0)  THEN
     999      IF (N == 2)  THEN
     1000        ! NO CONDITION AT LEFT END AND N = 2.
     1001        WK(2, 1) = ONE
     1002        WK(1, 1) = ONE
     1003        D(1, 1) = TWO * WK(2, 2)
     1004      ELSE
     1005        ! NOT-A-KNOT CONDITION AT LEFT END AND N .GT. 2.
     1006        WK(2, 1) = WK(1, 3)
     1007        WK(1, 1) = WK(1, 2) + WK(1, 3)
     1008        D(1, 1) = ((WK(1, 2) + TWO * WK(1, 1)) * WK(2, 2) * WK(1, 3) &
     1009                + WK(1, 2)**2 * WK(2, 3)) / WK(1, 1)
     1010      ENDIF
     1011    ELSE IF (IBEG == 1)  THEN
     1012      ! SLOPE PRESCRIBED AT LEFT END.
     1013      WK(2, 1) = ONE
     1014      WK(1, 1) = ZERO
     1015    ELSE
     1016      ! SECOND DERIVATIVE PRESCRIBED AT LEFT END.
     1017      WK(2, 1) = TWO
     1018      WK(1, 1) = ONE
     1019      D(1, 1) = THREE * WK(2, 2) - HALF * WK(1, 2) * D(1, 1)
     1020    ENDIF
     1021    !
     1022    !  IF THERE ARE INTERIOR KNOTS, GENERATE THE CORRESPONDING EQUATIONS AND
     1023    !  CARRY OUT THE FORWARD PASS OF GAUSS ELIMINATION, AFTER WHICH THE J-TH
     1024    !  EQUATION READS    WK(2,J)*S(J) + WK(1,J)*S(J+1) = D(1,J).
     1025    !
     1026    NM1 = N - 1
     1027    IF (NM1 > 1)  THEN
     1028      DO J = 2, NM1
     1029        IF (WK(2, J - 1) == ZERO)  GO TO 5008
     1030        G = -WK(1, J + 1) / WK(2, J - 1)
     1031        D(1, J) = G * D(1, J - 1) &
     1032                + THREE * (WK(1, J) * WK(2, J + 1) + WK(1, J + 1) * WK(2, J))
     1033        WK(2, J) = G * WK(1, J - 1) + TWO * (WK(1, J) + WK(1, J + 1))
     1034      END DO
     1035    ENDIF
     1036    !
     1037    !  CONSTRUCT LAST EQUATION FROM SECOND BOUNDARY CONDITION, OF THE FORM
     1038    !       (-G*WK(2,N-1))*S(N-1) + WK(2,N)*S(N) = D(1,N)
     1039    !
     1040    ! IF SLOPE IS PRESCRIBED AT RIGHT END, ONE CAN GO DIRECTLY TO BACK-
     1041    ! SUBSTITUTION, SINCE ARRAYS HAPPEN TO BE SET UP JUST RIGHT FOR IT
     1042    ! AT THIS POINT.
     1043    IF (IEND == 1)  GO TO 30
     1044    !
     1045    IF (IEND == 0)  THEN
     1046      IF (N==2 .AND. IBEG==0)  THEN
     1047        ! NOT-A-KNOT AT RIGHT ENDPOINT AND AT LEFT ENDPOINT AND N = 2.
     1048        D(1, 2) = WK(2, 2)
     1049        GO TO 30
     1050      ELSE IF ((N==2) .OR. (N==3 .AND. IBEG==0))  THEN
     1051        ! EITHER (N=3 AND NOT-A-KNOT ALSO AT LEFT) OR (N=2 AND *NOT*
     1052        ! NOT-A-KNOT AT LEFT END POINT).
     1053        D(1, N) = TWO * WK(2, N)
     1054        WK(2, N) = ONE
     1055        IF (WK(2, N - 1) == ZERO)  GO TO 5008
     1056        G = -ONE / WK(2, N - 1)
     1057      ELSE
     1058        ! NOT-A-KNOT AND N .GE. 3, AND EITHER N.GT.3 OR  ALSO NOT-A-
     1059        ! KNOT AT LEFT END POINT.
     1060        G = WK(1, N - 1) + WK(1, N)
     1061        ! DO NOT NEED TO CHECK FOLLOWING DENOMINATORS (X-DIFFERENCES).
     1062        D(1, N) = ((WK(1, N) + TWO * G) * WK(2, N) * WK(1, N - 1) &
     1063                + WK(1, N)**2 * (F(1, N - 1) - F(1, N - 2)) / WK(1, N - 1)) / G
     1064        IF (WK(2, N - 1) == ZERO)  GO TO 5008
     1065        G = -G / WK(2, N - 1)
     1066        WK(2, N) = WK(1, N - 1)
     1067      ENDIF
     1068    ELSE
     1069      ! SECOND DERIVATIVE PRESCRIBED AT RIGHT ENDPOINT.
     1070      D(1, N) = THREE * WK(2, N) + HALF * WK(1, N) * D(1, N)
     1071      WK(2, N) = TWO
     1072      IF (WK(2, N - 1) == ZERO)  GO TO 5008
     1073      G = -ONE / WK(2, N - 1)
     1074    ENDIF
     1075    !
     1076    !  COMPLETE FORWARD PASS OF GAUSS ELIMINATION.
     1077    !
     1078    WK(2, N) = G * WK(1, N - 1) + WK(2, N)
     1079    IF (WK(2, N) == ZERO)   GO TO 5008
     1080    D(1, N) = (G * D(1, N - 1) + D(1, N)) / WK(2, N)
     1081    !
     1082    !  CARRY OUT BACK SUBSTITUTION
     1083    !
     1084    30   CONTINUE
     1085    DO J = NM1, 1, -1
     1086      IF (WK(2, J) == ZERO)  GO TO 5008
     1087      D(1, J) = (D(1, J) - WK(1, J) * D(1, J + 1)) / WK(2, J)
     1088    END DO
     1089    ! --------------------(  END  CODING FROM CUBSPL )--------------------
     1090    !
     1091    !  NORMAL RETURN.
     1092    !
     1093    RETURN
     1094    !
     1095    !  ERROR RETURNS.
     1096    !
     1097    5001   CONTINUE
     1098    ! N.LT.2 RETURN.
     1099    IERR = -1
     1100    CALL XERMSG ('SLATEC', 'PCHSP', &
     1101            'NUMBER OF DATA POINTS LESS THAN TWO', IERR, 1)
     1102    RETURN
     1103    !
     1104    5002   CONTINUE
     1105    ! INCFD.LT.1 RETURN.
     1106    IERR = -2
     1107    CALL XERMSG ('SLATEC', 'PCHSP', 'INCREMENT LESS THAN ONE', IERR, &
     1108            1)
     1109    RETURN
     1110    !
     1111    5003   CONTINUE
     1112    ! X-ARRAY NOT STRICTLY INCREASING.
     1113    IERR = -3
     1114    CALL XERMSG ('SLATEC', 'PCHSP', 'X-ARRAY NOT STRICTLY INCREASING' &
     1115            , IERR, 1)
     1116    RETURN
     1117    !
     1118    5004   CONTINUE
     1119    ! IC OUT OF RANGE RETURN.
     1120    IERR = IERR - 3
     1121    CALL XERMSG ('SLATEC', 'PCHSP', 'IC OUT OF RANGE', IERR, 1)
     1122    RETURN
     1123    !
     1124    5007   CONTINUE
     1125    ! NWK TOO SMALL RETURN.
     1126    IERR = -7
     1127    CALL XERMSG ('SLATEC', 'PCHSP', 'WORK ARRAY TOO SMALL', IERR, 1)
     1128    RETURN
     1129    !
     1130    5008   CONTINUE
     1131    ! SINGULAR SYSTEM.
     1132    !   *** THEORETICALLY, THIS CAN ONLY OCCUR IF SUCCESSIVE X-VALUES   ***
     1133    !   *** ARE EQUAL, WHICH SHOULD ALREADY HAVE BEEN CAUGHT (IERR=-3). ***
     1134    IERR = -8
     1135    CALL XERMSG ('SLATEC', 'PCHSP', 'SINGULAR LINEAR SYSTEM', IERR, &
     1136            1)
     1137    RETURN
     1138    !
     1139    5009   CONTINUE
     1140    ! ERROR RETURN FROM PCHDF.
     1141    !   *** THIS CASE SHOULD NEVER OCCUR ***
     1142    IERR = -9
     1143    CALL XERMSG ('SLATEC', 'PCHSP', 'ERROR RETURN FROM PCHDF', IERR, &
     1144            1)
     1145    RETURN
     1146    !------------- LAST LINE OF PCHSP FOLLOWS ------------------------------
     1147  END SUBROUTINE PCHSP
     1148
     1149
     1150END MODULE lmdz_libmath_pch

  • LMDZ6/branches/Amaury_dev/libf/misc/new_unit_m.F90

    r5103 r5113  
    11module new_unit_m
    22
    3   implicit none
     3  IMPLICIT NONE
    44
    55contains

  • LMDZ6/branches/Amaury_dev/libf/misc/nrtype.F90

    r2228 r5113  
    11MODULE nrtype
    22
    3   implicit none
     3  IMPLICIT NONE
    44
    55  integer, parameter:: k8 = selected_real_kind(13)

  • LMDZ6/branches/Amaury_dev/libf/misc/regr1_step_av_m.F90

    r5101 r5113  
    44  ! Author: Lionel GUEZ
    55
    6   implicit none
     6  IMPLICIT NONE
    77
    88  interface regr1_step_av
    … …  
    2121  end interface
    2222
    23   private
     23  PRIVATE
    2424  public regr1_step_av
    2525
    … …  
    3030    ! "vs" has rank 1.
    3131
    32     use assert_eq_m, only: assert_eq
    33     use assert_m, only: assert
     32    use lmdz_assert_eq, only: assert_eq
     33    use lmdz_assert, only: assert
    3434    use interpolation, only: locate
    3535
    … …  
    8989    ! "vs" has rank 2.
    9090
    91     use assert_eq_m, only: assert_eq
    92     use assert_m, only: assert
     91    use lmdz_assert_eq, only: assert_eq
     92    use lmdz_assert, only: assert
    9393    use interpolation, only: locate
    9494
    … …  
    149149    ! "vs" has rank 3.
    150150
    151     use assert_eq_m, only: assert_eq
    152     use assert_m, only: assert
     151    use lmdz_assert_eq, only: assert_eq
     152    use lmdz_assert, only: assert
    153153    use interpolation, only: locate
    154154
    … …  
    210210    ! "vs" has rank 4.
    211211
    212     use assert_eq_m, only: assert_eq
    213     use assert_m, only: assert
     212    use lmdz_assert_eq, only: assert_eq
     213    use lmdz_assert, only: assert
    214214    use interpolation, only: locate
    215215

  • LMDZ6/branches/Amaury_dev/libf/misc/regr_conserv_m.F90

    r5099 r5113  
    11MODULE regr_conserv_m
    22
    3   USE assert_eq_m,   ONLY: assert_eq
    4   USE assert_m,      ONLY: assert
     3  USE lmdz_assert_eq,   ONLY: assert_eq
     4  USE lmdz_assert,      ONLY: assert
    55  USE interpolation, ONLY: locate
    66

  • LMDZ6/branches/Amaury_dev/libf/misc/regr_lint_m.F90

    r5101 r5113  
    11MODULE regr_lint_m
    22
    3   USE assert_eq_m,   ONLY: assert_eq
    4   USE assert_m,      ONLY: assert
     3  USE lmdz_assert_eq,   ONLY: assert_eq
     4  USE lmdz_assert,      ONLY: assert
    55  USE interpolation, ONLY: hunt
    66

  • LMDZ6/branches/Amaury_dev/libf/misc/vampir.F90

    r5103 r5113  
    1616
    1717  SUBROUTINE InitVampir
    18     implicit none
     18    IMPLICIT NONE
    1919
    2020#ifdef USE_VT
    … …  
    5050
    5151  SUBROUTINE VTb(number)
    52     implicit none
     52    IMPLICIT NONE
    5353    INTEGER :: number
    5454#ifdef USE_VT   
    … …  
    6767
    6868  SUBROUTINE VTe(number)
    69     implicit none
     69    IMPLICIT NONE
    7070    INTEGER :: Number
    7171#ifdef USE_VT   

  • LMDZ6/branches/Amaury_dev/libf/misc/write_field.F90

    r5103 r5113  
    2222
    2323  function GetFieldIndex(name)
    24     implicit none
     24    IMPLICIT NONE
    2525    integer :: GetFieldindex
    2626    character(len = *) :: name
    … …  
    4141
    4242  subroutine WriteField3d(name, Field)
    43     implicit none
     43    IMPLICIT NONE
    4444    character(len = *) :: name
    4545    real, dimension(:, :, :) :: Field
    … …  
    5252
    5353  subroutine WriteField2d(name, Field)
    54     implicit none
     54    IMPLICIT NONE
    5555    character(len = *) :: name
    5656    real, dimension(:, :) :: Field
    … …  
    6363
    6464  subroutine WriteField1d(name, Field)
    65     implicit none
     65    IMPLICIT NONE
    6666    character(len = *) :: name
    6767    real, dimension(:) :: Field
    … …  
    7474
    7575  subroutine WriteField_gen(name, Field, dimx, dimy, dimz)
    76     implicit none
     76    IMPLICIT NONE
    7777    character(len = *) :: name
    7878    integer :: dimx, dimy, dimz
    … …  
    108108
    109109  subroutine CreateNewField(name, dimx, dimy, dimz)
    110     implicit none
     110    IMPLICIT NONE
    111111    character(len = *) :: name
    112112    integer :: dimx, dimy, dimz
    … …  
    129129
    130130  subroutine write_field1D(name, Field)
    131     implicit none
     131    IMPLICIT NONE
    132132
    133133    integer, parameter :: MaxDim = 1
    … …  
    170170
    171171  subroutine write_field2D(name, Field)
    172     implicit none
     172    IMPLICIT NONE
    173173
    174174    integer, parameter :: MaxDim = 2
    … …  
    229229
    230230  subroutine write_field3D(name, Field)
    231     implicit none
     231    IMPLICIT NONE
    232232
    233233    integer, parameter :: MaxDim = 3
Note: See TracChangeset for help on using the changeset viewer.