1 | MODULE CPLEDN_MOD |
---|
2 | CONTAINS |
---|
3 | SUBROUTINE CPLEDN(KN,KDBLE,PX,DDX,KFLAG,PW,PXN,DDXN,PXMOD) |
---|
4 | |
---|
5 | !**** *CPLEDN* - Routine to compute the Legendre polynomial of degree N |
---|
6 | |
---|
7 | ! Purpose. |
---|
8 | ! -------- |
---|
9 | ! Computes Legendre polynomial of degree N |
---|
10 | |
---|
11 | !** Interface. |
---|
12 | ! ---------- |
---|
13 | ! *CALL* *CPLEDN(KN,KDBLE,PX,DDX,KFLAG,PW,PXN,DDXN,PXMOD)* |
---|
14 | |
---|
15 | ! Explicit arguments : |
---|
16 | ! -------------------- |
---|
17 | ! KN : Degree of the Legendre polynomial |
---|
18 | ! KDBLE : 0, single precision |
---|
19 | ! 1, double precision |
---|
20 | ! PX : abcissa where the computations are performed |
---|
21 | ! DDX : id in double precision |
---|
22 | ! KFLAG : When KFLAG.EQ.1 computes the weights |
---|
23 | ! PW : Weight of the quadrature at PXN |
---|
24 | ! PXN : new abscissa (Newton iteration) |
---|
25 | ! DDXN : id in double precision |
---|
26 | ! PXMOD : PXN-PX |
---|
27 | |
---|
28 | ! Implicit arguments : |
---|
29 | ! -------------------- |
---|
30 | ! None |
---|
31 | |
---|
32 | ! Method. |
---|
33 | ! ------- |
---|
34 | ! See documentation |
---|
35 | |
---|
36 | ! Externals. |
---|
37 | ! ---------- |
---|
38 | ! None |
---|
39 | |
---|
40 | ! Reference. |
---|
41 | ! ---------- |
---|
42 | ! ECMWF Research Department documentation of the IFS |
---|
43 | |
---|
44 | ! Author. |
---|
45 | ! ------- |
---|
46 | ! Mats Hamrud and Philippe Courtier *ECMWF* |
---|
47 | |
---|
48 | ! Modifications. |
---|
49 | ! -------------- |
---|
50 | ! Original : 87-10-15 |
---|
51 | ! Michel Rochas, 90-08-30 (Lobatto+cleaning) |
---|
52 | ! ------------------------------------------------------------------ |
---|
53 | |
---|
54 | |
---|
55 | |
---|
56 | USE PARKIND1 ,ONLY : JPIM ,JPRB |
---|
57 | USE PARKIND2 ,ONLY : JPRH |
---|
58 | |
---|
59 | IMPLICIT NONE |
---|
60 | |
---|
61 | |
---|
62 | ! DUMMY INTEGER SCALARS |
---|
63 | INTEGER(KIND=JPIM) :: KDBLE |
---|
64 | INTEGER(KIND=JPIM) :: KFLAG |
---|
65 | INTEGER(KIND=JPIM) :: KN |
---|
66 | |
---|
67 | ! DUMMY REAL SCALARS |
---|
68 | REAL(KIND=JPRB) :: PW |
---|
69 | REAL(KIND=JPRB) :: PX |
---|
70 | REAL(KIND=JPRB) :: PXMOD |
---|
71 | REAL(KIND=JPRB) :: PXN |
---|
72 | |
---|
73 | |
---|
74 | REAL(KIND=JPRH) :: DDX,DDXN,DLX,DLK,DLKM1,DLKM2,DLLDN,DLXN,DLMOD |
---|
75 | REAL(KIND=JPRH) :: DLG,DLGDN |
---|
76 | |
---|
77 | INTEGER(KIND=JPIM), PARAMETER :: JPKS=KIND(PX) |
---|
78 | INTEGER(KIND=JPIM), PARAMETER :: JPKD=KIND(DDX) |
---|
79 | |
---|
80 | ! LOCAL INTEGER SCALARS |
---|
81 | INTEGER(KIND=JPIM) :: IZN, JN |
---|
82 | |
---|
83 | ! LOCAL REAL SCALARS |
---|
84 | REAL(KIND=JPRB) :: ZG, ZGDN, ZK, ZKM1, ZKM2, ZLDN, ZMOD, ZX, ZXN |
---|
85 | |
---|
86 | |
---|
87 | ! ----------------------------------------------------------------- |
---|
88 | |
---|
89 | !* 1. Single precision computations. |
---|
90 | ! ------------------------------ |
---|
91 | |
---|
92 | IZN = KN |
---|
93 | |
---|
94 | ZK = 0.0_JPRB |
---|
95 | DLK = 0.0_JPRB |
---|
96 | DLXN = 0.0_JPRB |
---|
97 | IF(KDBLE == 0)THEN |
---|
98 | |
---|
99 | !* 1.1 NEWTON ITERATION STEP. |
---|
100 | |
---|
101 | ZKM2 = 1 |
---|
102 | ZKM1 = PX |
---|
103 | ZX = PX |
---|
104 | DO JN=2,IZN |
---|
105 | ZK = (REAL(2*JN-1,JPRB)*ZX*ZKM1-REAL(JN-1,JPRB)*ZKM2)/REAL(JN,JPRB) |
---|
106 | ZKM2 = ZKM1 |
---|
107 | ZKM1 = ZK |
---|
108 | ENDDO |
---|
109 | ZKM1 = ZKM2 |
---|
110 | ZLDN = (REAL(KN,JPRB)*(ZKM1-ZX*ZK))/(1.0_JPRB-ZX*ZX) |
---|
111 | ZMOD = -ZK/ZLDN |
---|
112 | ZXN = ZX+ZMOD |
---|
113 | PXN = ZXN |
---|
114 | DDXN = REAL(ZXN,JPKD) |
---|
115 | PXMOD = ZMOD |
---|
116 | |
---|
117 | ! ------------------------------------------------------------------ |
---|
118 | |
---|
119 | !* 2. Double precision computations. |
---|
120 | ! ------------------------------ |
---|
121 | |
---|
122 | ELSE |
---|
123 | |
---|
124 | !* 2.1 NEWTON ITERATION STEP. |
---|
125 | |
---|
126 | DLKM2 = 1.0_JPRB |
---|
127 | DLKM1 = DDX |
---|
128 | DLX = DDX |
---|
129 | DO JN=2,IZN |
---|
130 | DLK = (REAL(2*JN-1,JPKD)*DLX*DLKM1-REAL(JN-1,JPKD)*DLKM2)/REAL(JN,JPKD) |
---|
131 | DLKM2 = DLKM1 |
---|
132 | DLKM1 = DLK |
---|
133 | ENDDO |
---|
134 | DLKM1 = DLKM2 |
---|
135 | DLLDN = (REAL(KN,JPKD)*(DLKM1-DLX*DLK))/(1.0_JPRB-DLX*DLX) |
---|
136 | DLMOD = -DLK/DLLDN |
---|
137 | DLXN = DLX+DLMOD |
---|
138 | PXN = REAL(DLXN,JPKS) |
---|
139 | DDXN = DLXN |
---|
140 | PXMOD = REAL(DLMOD,JPKS) |
---|
141 | ENDIF |
---|
142 | ! ------------------------------------------------------------------ |
---|
143 | |
---|
144 | !* 3. Computes weight. |
---|
145 | ! ---------------- |
---|
146 | |
---|
147 | |
---|
148 | IF(KFLAG == 1)THEN |
---|
149 | DLKM2 = 1.0_JPRB |
---|
150 | DLKM1 = DLXN |
---|
151 | DLX = DLXN |
---|
152 | DO JN=2,IZN |
---|
153 | DLK = (REAL(2*JN-1,JPKD)*DLX*DLKM1-REAL(JN-1,JPKD)*DLKM2)/REAL(JN,JPKD) |
---|
154 | DLKM2 = DLKM1 |
---|
155 | DLKM1 = DLK |
---|
156 | ENDDO |
---|
157 | DLKM1 = DLKM2 |
---|
158 | PW = REAL((1.0_JPRB-DLX*DLX)/(REAL(KN*KN,JPKD)*DLKM1*DLKM1),JPKS) |
---|
159 | ENDIF |
---|
160 | |
---|
161 | ! ------------------------------------------------------------------ |
---|
162 | |
---|
163 | END SUBROUTINE CPLEDN |
---|
164 | END MODULE CPLEDN_MOD |
---|