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 |