*IF DEF,C93_1A,OR,DEF,C93_2A GNF0F402.10
C ******************************COPYRIGHT****************************** GTS2F400.2557
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved. GTS2F400.2558
C GTS2F400.2559
C Use, duplication or disclosure of this code is subject to the GTS2F400.2560
C restrictions as set forth in the contract. GTS2F400.2561
C GTS2F400.2562
C Meteorological Office GTS2F400.2563
C London Road GTS2F400.2564
C BRACKNELL GTS2F400.2565
C Berkshire UK GTS2F400.2566
C RG12 2SZ GTS2F400.2567
C GTS2F400.2568
C If no contract has been raised with this copy of the code, the use, GTS2F400.2569
C duplication or disclosure of it is strictly prohibited. Permission GTS2F400.2570
C to do so must first be obtained in writing from the Head of Numerical GTS2F400.2571
C Modelling at the above address. GTS2F400.2572
C ******************************COPYRIGHT****************************** GTS2F400.2573
C GTS2F400.2574
CLL SUBROUTINE EVAL_SP------------------------------------------------- EVALSP1A.3
CLL EVALSP1A.4
CLL PURPOSE: To evaluate a cubic spline at a given point. EVALSP1A.5
CLL Outputs spline value,1st and 2nd derivatives. EVALSP1A.6
CLL Can be called after subroutine SPLINE which outputs the EVALSP1A.7
CLL second derivative at all data points. EVALSP1A.8
CLL Tested under compiler CFT77 EVALSP1A.9
CLL Tested under OS version 5.1 EVALSP1A.10
CLL EVALSP1A.11
CLL Author J.T.Heming EVALSP1A.12
CLL EVALSP1A.13
CLL Model Modification history from model version 3.0: EVALSP1A.14
CLL version Date EVALSP1A.15
CLL vn4.2 07/11/96:Allow this routine to be used in C93_2A (Farnon) GNF0F402.9
CLL EVALSP1A.16
CLL Logical components covered :D413 EVALSP1A.17
CLL EVALSP1A.18
CLL Project TASK: EVALSP1A.19
CLL EVALSP1A.20
CLL Programming standard: U M DOC Paper NO. 4, EVALSP1A.21
CLL EVALSP1A.22
CLL External documentation EVALSP1A.23
CLL EVALSP1A.24
CLLEND------------------------------------------------------------------ EVALSP1A.25
C EVALSP1A.26
C*L ARGUMENTS:--------------------------------------------------------- EVALSP1A.27
SUBROUTINE EVAL_SP( 3EVALSP1A.28
C data in EVALSP1A.29
& X,Y,AMU,H_POINTS,V_POINTS,X0, EVALSP1A.30
C data out EVALSP1A.31
& Y0,DERIV1,DERIV2) EVALSP1A.32
C*L EVALSP1A.33
C----------------------------------------------------------------------- EVALSP1A.34
IMPLICIT NONE EVALSP1A.35
C*---------------------------------------------------------------------- EVALSP1A.36
INTEGER EVALSP1A.37
* H_POINTS ! The number of splines to be evaluated EVALSP1A.38
*,V_POINTS ! The number of points in each spline profile EVALSP1A.39
C----------------------------------------------------------------------- EVALSP1A.40
REAL EVALSP1A.41
* X(V_POINTS) ! IN Independent variable array input data EVALSP1A.42
*,Y(H_POINTS,V_POINTS) ! IN Dependent variable array input data EVALSP1A.43
*,AMU(H_POINTS,V_POINTS) ! IN Array of 2nd derivatives EVALSP1A.44
*,X0 ! IN X-value at which spline is evaluated EVALSP1A.45
*,Y0(H_POINTS) ! OUT Array of spline values EVALSP1A.46
*,DERIV1(H_POINTS) ! OUT Array of 1st derivatives EVALSP1A.47
*,DERIV2(H_POINTS) ! OUT Array of 2nd derivatives EVALSP1A.48
C* EVALSP1A.49
C*L EVALSP1A.50
C----------------------------------------------------------------------- EVALSP1A.51
C Local Variables EVALSP1A.52
C----------------------------------------------------------------------- EVALSP1A.53
INTEGER EVALSP1A.54
* I,K,J ! Loop counters EVALSP1A.55
C----------------------------------------------------------------------- EVALSP1A.56
REAL EVALSP1A.57
* XD1,XD2,XD3,RXD3,XD12,RXD3O6 ! Work variables EVALSP1A.58
*,WORK1(H_POINTS) ! Work array EVALSP1A.59
C----------------------------------------------------------------------- EVALSP1A.60
LOGICAL EVALSP1A.61
* FOUND ! Set true once spline evaluated EVALSP1A.62
C----------------------------------------------------------------------- EVALSP1A.63
CL Check that X0 lies within range of X-values EVALSP1A.64
CL N.B. This routine assumes X decreases with height EVALSP1A.65
C----------------------------------------------------------------------- EVALSP1A.66
IF(X0.LE.X(1).AND.X0.GE.X(V_POINTS))THEN EVALSP1A.67
C----------------------------------------------------------------------- EVALSP1A.68
FOUND=.FALSE. EVALSP1A.69
C----------------------------------------------------------------------- EVALSP1A.70
CL Find adjacent X values between which X0 lies EVALSP1A.71
C----------------------------------------------------------------------- EVALSP1A.72
DO J=1,V_POINTS-1 EVALSP1A.73
IF(X0.GT.X(J+1).AND.(.NOT.FOUND))THEN EVALSP1A.74
FOUND=.TRUE. EVALSP1A.75
K=J EVALSP1A.76
ELSEIF(X0.EQ.X(V_POINTS).AND.(.NOT.FOUND))THEN EVALSP1A.77
FOUND=.TRUE. EVALSP1A.78
K=V_POINTS-1 EVALSP1A.79
ENDIF EVALSP1A.80
ENDDO EVALSP1A.81
C----------------------------------------------------------------------- EVALSP1A.82
CL Set local variables EVALSP1A.83
C----------------------------------------------------------------------- EVALSP1A.84
XD1=X(K+1)-X0 EVALSP1A.85
XD2=X0-X(K) EVALSP1A.86
XD3=X(K+1)-X(K) EVALSP1A.87
RXD3=1.0/XD3 EVALSP1A.88
RXD3O6=RXD3/6.0 EVALSP1A.89
XD12=XD1*XD2/6.0 EVALSP1A.90
C----------------------------------------------------------------------- EVALSP1A.91
CL Loop through horizontal field EVALSP1A.92
C----------------------------------------------------------------------- EVALSP1A.93
DO I=1,H_POINTS EVALSP1A.94
C----------------------------------------------------------------------- EVALSP1A.95
CL Calculate factor used in calculation of Y0 and DERIV1 EVALSP1A.96
C----------------------------------------------------------------------- EVALSP1A.97
WORK1(I)=(AMU(I,K+1)*(XD2+XD3)+AMU(I,K)*(XD1+XD3)) EVALSP1A.98
* *RXD3O6 EVALSP1A.99
C----------------------------------------------------------------------- EVALSP1A.100
CL Evaluate spline value Y0 at point X0 EVALSP1A.101
C----------------------------------------------------------------------- EVALSP1A.102
Y0(I)=((Y(I,K+1)*XD2)+(Y(I,K)*XD1))*RXD3 EVALSP1A.103
Y0(I)=Y0(I)-WORK1(I)*XD1*XD2 EVALSP1A.104
C----------------------------------------------------------------------- EVALSP1A.105
CL Evaluate first derivative at point X0 EVALSP1A.106
C----------------------------------------------------------------------- EVALSP1A.107
DERIV1(I)=(Y(I,K+1)-Y(I,K)+(AMU(I,K+1)-AMU(I,K))*XD12) EVALSP1A.108
* *RXD3 EVALSP1A.109
DERIV1(I)=DERIV1(I)+WORK1(I)*(XD2-XD1) EVALSP1A.110
C----------------------------------------------------------------------- EVALSP1A.111
CL Evaluate second deriavtive at point X0 EVALSP1A.112
C----------------------------------------------------------------------- EVALSP1A.113
DERIV2(I)=((AMU(I,K+1)*XD1)+(AMU(I,K)*XD2))*RXD3 EVALSP1A.114
ENDDO EVALSP1A.115
C----------------------------------------------------------------------- EVALSP1A.116
C Error message EVALSP1A.117
C----------------------------------------------------------------------- EVALSP1A.118
ELSE EVALSP1A.119
WRITE(6,999)X0 EVALSP1A.120
999 FORMAT(' Spline not evaluated - X0 out of range. X0=',F10.4) EVALSP1A.121
ENDIF EVALSP1A.122
C======================================================================= EVALSP1A.123
C END OF SUBROUTINE EVAL_SP EVALSP1A.124
C======================================================================= EVALSP1A.125
RETURN EVALSP1A.126
END EVALSP1A.127
C======================================================================= EVALSP1A.128
*ENDIF EVALSP1A.129