1 |
|
|
MODULE SULEG_MOD |
2 |
|
|
CONTAINS |
3 |
|
|
SUBROUTINE SULEG |
4 |
|
|
|
5 |
|
|
USE PARKIND1 ,ONLY : JPIM ,JPRB |
6 |
|
|
USE PARKIND2 ,ONLY : JPRH |
7 |
|
|
|
8 |
|
|
USE TPM_GEN |
9 |
|
|
USE TPM_DIM |
10 |
|
|
USE TPM_CONSTANTS |
11 |
|
|
USE TPM_DISTR |
12 |
|
|
USE TPM_FIELDS |
13 |
|
|
|
14 |
|
|
!USE SUGAW_MOD |
15 |
|
|
USE SUPOL_MOD |
16 |
|
|
USE SUTRLE_MOD |
17 |
|
|
|
18 |
|
|
#ifdef DOC |
19 |
|
|
|
20 |
|
|
!**** *SULEG * - initialize the Legendre polynomials |
21 |
|
|
|
22 |
|
|
! Purpose. |
23 |
|
|
! -------- |
24 |
|
|
! Initialize COMMON YOMLEG |
25 |
|
|
|
26 |
|
|
!** Interface. |
27 |
|
|
! ---------- |
28 |
|
|
! *CALL* *SULEG* |
29 |
|
|
|
30 |
|
|
! Explicit arguments : |
31 |
|
|
! -------------------- |
32 |
|
|
|
33 |
|
|
! Implicit arguments : |
34 |
|
|
! -------------------- |
35 |
|
|
! COMMON YOMLEG |
36 |
|
|
|
37 |
|
|
! Method. |
38 |
|
|
! ------- |
39 |
|
|
! See documentation |
40 |
|
|
|
41 |
|
|
! Externals. |
42 |
|
|
! ---------- |
43 |
|
|
! SUGAW (Gaussian latitudes) |
44 |
|
|
! SUPOLM (polynomials) |
45 |
|
|
! LFI routines for external IO's |
46 |
|
|
! Called by SUGEM. |
47 |
|
|
|
48 |
|
|
! Reference. |
49 |
|
|
! ---------- |
50 |
|
|
! ECMWF Research Department documentation of the IFS |
51 |
|
|
|
52 |
|
|
! Author. |
53 |
|
|
! ------- |
54 |
|
|
! Mats Hamrud and Philippe Courtier *ECMWF* |
55 |
|
|
|
56 |
|
|
! Modifications. |
57 |
|
|
! -------------- |
58 |
|
|
! Original : 87-10-15 |
59 |
|
|
! MODIFICATION : 91-04 J.M. Piriou: |
60 |
|
|
! - Read gaussian latitudes and PNM on LFI |
61 |
|
|
! - If file missing, computes |
62 |
|
|
! 91-04 M.Hamrud: |
63 |
|
|
! - IO Scheme introduced |
64 |
|
|
! MODIFICATION : 91-07-03 P.Courtier suppress derivatives |
65 |
|
|
! MODIFICATION : 91-07-03 P.Courtier computes RATATH and RACTHE |
66 |
|
|
! MODIFICATION : 91-07-03 P.Courtier change upper limit (NSMAX+1) |
67 |
|
|
! MODIFICATION : 91-07-03 P.Courtier change ordering |
68 |
|
|
! Order of the PNM in the file, as in the model : |
69 |
|
|
! - increasing wave numbers m |
70 |
|
|
! - for a given m, from n=NSMAX+1 to m |
71 |
|
|
! MODIFICATION : 92-07-02 R. Bubnova: shift RATATH calculation |
72 |
|
|
! to SUGEM1 |
73 |
|
|
! MODIFICATION : 92-12-17 P.Courtier multitask computations |
74 |
|
|
! Modified by R. EL Khatib : 93-04-02 Set-up defaults controled by LECMWF |
75 |
|
|
! MODIFICATION : 93-03-19 D.Giard : n <= NTMAX |
76 |
|
|
! K. YESSAD : 93-05-11 : DLMU --> global array DRMU(NDGSA:NDGEN). |
77 |
|
|
! (not stored currently on LFI files). |
78 |
|
|
! MODIFICATION : 94-02-03 R. El Khatib : subroutine SULEG2 to write out |
79 |
|
|
! the Leg. polynomials on workfile or LFI file |
80 |
|
|
! Modification : 94-08-31 M. Tolstykh: Setup for CUD interpolation |
81 |
|
|
! Modified by K. YESSAD (MARCH 1995): Extra-latitudes computations |
82 |
|
|
! according to value of NDGSUR and LRPOLE only. |
83 |
|
|
! + change fancy loop numbering. |
84 |
|
|
! Modified 98-08-10 by K. YESSAD: removal of LRPOLE option. |
85 |
|
|
! - removal of LRPOLE in YOMCT0. |
86 |
|
|
! - removal of code under LRPOLE. |
87 |
|
|
! ------------------------------------------------------------------ |
88 |
|
|
|
89 |
|
|
#endif |
90 |
|
|
|
91 |
|
|
|
92 |
|
|
IMPLICIT NONE |
93 |
|
|
|
94 |
|
|
|
95 |
|
|
! ------------------------------------------------------------------ |
96 |
|
|
REAL(KIND=JPRB),ALLOCATABLE :: ZPNMG(:,:) |
97 |
|
|
|
98 |
|
|
REAL(KIND=JPRH) :: DLRMU(R%NDGL) |
99 |
|
|
REAL(KIND=JPRH) :: DLC(0:R%NTMAX+1,0:R%NTMAX+1) |
100 |
|
|
REAL(KIND=JPRH) :: DLD(0:R%NTMAX+1,0:R%NTMAX+1) |
101 |
|
|
REAL(KIND=JPRH) :: DLE(0:R%NTMAX+1,0:R%NTMAX+1) |
102 |
|
|
REAL(KIND=JPRH) :: DLA(0:R%NTMAX+1),DLB(0:R%NTMAX+1),DLF(0:R%NTMAX+1) |
103 |
|
|
REAL(KIND=JPRH) :: DLG(0:R%NTMAX+1),DLH(0:R%NTMAX+1),DLI(0:R%NTMAX+1) |
104 |
|
|
REAL(KIND=JPRH) :: DLPOL(0:R%NTMAX+1,0:R%NTMAX+1) |
105 |
|
|
! ------------------------------------------------------------------ |
106 |
|
|
|
107 |
|
|
INTEGER(KIND=JPIM), PARAMETER :: JPKS=KIND(ZPNMG) |
108 |
|
|
INTEGER(KIND=JPIM), PARAMETER :: JPKD=KIND(DLG) |
109 |
|
|
|
110 |
|
|
! ------------------------------------------------------------------ |
111 |
|
|
REAL(KIND=JPRH) :: DA,DC,DD,DE |
112 |
|
|
INTEGER(KIND=JPIM) :: KKN, KKM |
113 |
|
|
|
114 |
|
|
! LOCAL |
115 |
|
|
INTEGER(KIND=JPIM) :: IGLLOC, INM, IM , ICOUNT,& |
116 |
|
|
&JGL, JM, JMLOC, JN, JNM |
117 |
|
|
|
118 |
|
|
|
119 |
|
|
LOGICAL :: LLP1,LLP2 |
120 |
|
|
|
121 |
|
|
|
122 |
|
|
DC(KKN,KKM)=SQRT( (REAL(2*KKN+1,JPKD)*REAL(KKN+KKM-1,JPKD)& |
123 |
|
|
&*REAL(KKN+KKM-3,JPKD))& |
124 |
|
|
&/ (REAL(2*KKN-3,JPKD)*REAL(KKN+KKM,JPKD)& |
125 |
|
|
&*REAL(KKN+KKM-2,JPKD)) ) |
126 |
|
|
DD(KKN,KKM)=SQRT( (REAL(2*KKN+1,JPKD)*REAL(KKN+KKM-1,JPKD)& |
127 |
|
|
&*REAL(KKN-KKM+1,JPKD))& |
128 |
|
|
&/ (REAL(2*KKN-1,JPKD)*REAL(KKN+KKM,JPKD)& |
129 |
|
|
&*REAL(KKN+KKM-2,JPKD)) ) |
130 |
|
|
DE(KKN,KKM)=SQRT( (REAL(2*KKN+1,JPKD)*REAL(KKN-KKM,JPKD))& |
131 |
|
|
&/ (REAL(2*KKN-1,JPKD)*REAL(KKN+KKM,JPKD)) ) |
132 |
|
|
DA(KKN,KKM)=SQRT( (REAL(2*KKN+1,JPKD)*REAL(KKN-KKM,JPKD)& |
133 |
|
|
&*REAL(KKN+KKM,JPKD))& |
134 |
|
|
&/ REAL(2*KKN-1,JPKD) ) |
135 |
|
|
|
136 |
|
|
! ------------------------------------------------------------------ |
137 |
|
|
ALLOCATE(ZPNMG(R%NSPOLEG,D%NLEI3D)) |
138 |
|
|
|
139 |
|
|
!* 0. Some initializations. |
140 |
|
|
! --------------------- |
141 |
|
|
|
142 |
|
|
LLP1 = NPRINTLEV>0 |
143 |
|
|
LLP2 = NPRINTLEV>1 |
144 |
|
|
IF(LLP1) WRITE(NOUT,*) '=== ENTER ROUTINE SULEG ===' |
145 |
|
|
|
146 |
|
|
!CALL GSTATS(140,0) !MPL 4.12.08 |
147 |
|
|
ALLOCATE(F%RPNM(R%NLEI3,D%NSPOLEGL)) |
148 |
|
|
IF (LLP2) WRITE(NOUT,9) 'F%RPNM ',SIZE(F%RPNM),SHAPE(F%RPNM) |
149 |
|
|
ALLOCATE(F%RMU(R%NDGL)) |
150 |
|
|
IF (LLP2) WRITE(NOUT,9) 'F%RMU ',SIZE(F%RMU ),SHAPE(F%RMU ) |
151 |
|
|
ALLOCATE(F%RW(R%NDGL)) |
152 |
|
|
IF (LLP2) WRITE(NOUT,9) 'F%RW ',SIZE(F%RW ),SHAPE(F%RW ) |
153 |
|
|
ALLOCATE(F%R1MU2(R%NDGL)) |
154 |
|
|
IF (LLP2) WRITE(NOUT,9) 'F%R1MU2 ',SIZE(F%R1MU2),SHAPE(F%R1MU2 ) |
155 |
|
|
ALLOCATE(F%RACTHE(R%NDGL)) |
156 |
|
|
IF (LLP2) WRITE(NOUT,9) 'F%RACTHE ',SIZE(F%RACTHE),SHAPE(F%RACTHE ) |
157 |
|
|
|
158 |
|
|
!CALL GSTATS(1801,0) ! MPL 4.12.08 |
159 |
|
|
DO JNM=1,D%NSPOLEGL |
160 |
|
|
F%RPNM(R%NLEI3,JNM) = 0.0_JPRB |
161 |
|
|
ENDDO |
162 |
|
|
!CALL GSTATS(1801,1) ! MPL 4.12.08 |
163 |
|
|
|
164 |
|
|
! ------------------------------------------------------------------ |
165 |
|
|
|
166 |
|
|
!* 3.1 Gaussian latitudes and weights |
167 |
|
|
!CALL SUGAW(R%NDGL,F%RMU,DLRMU,F%RW) |
168 |
|
|
|
169 |
|
|
!* 3.2 Computes related arrays |
170 |
|
|
|
171 |
|
|
DO JGL=1,R%NDGL |
172 |
|
|
F%R1MU2(JGL) = REAL(1.0_JPRB-DLRMU(JGL)*DLRMU(JGL),JPKS) |
173 |
|
|
F%RACTHE(JGL) = REAL(1.0_JPRB/SQRT(1.0_JPRB-DLRMU(JGL)*DLRMU(JGL))/& |
174 |
|
|
&REAL(RA,JPKD),JPKS) |
175 |
|
|
ENDDO |
176 |
|
|
|
177 |
|
|
!* 3.3 Working arrays |
178 |
|
|
DO JN=3,R%NTMAX+1 |
179 |
|
|
DO JM=2,JN-1 |
180 |
|
|
DLC(JM,JN) = DC(JN,JM) |
181 |
|
|
DLD(JM,JN) = DD(JN,JM) |
182 |
|
|
DLE(JM,JN) = DE(JN,JM) |
183 |
|
|
ENDDO |
184 |
|
|
ENDDO |
185 |
|
|
|
186 |
|
|
DO JN=1,R%NTMAX+1 |
187 |
|
|
DLA(JN) = SQRT(REAL(2*JN+1,JPKD)) |
188 |
|
|
DLB(JN) = SQRT(REAL(2*JN+1,JPKD)/REAL(JN*(JN+1),JPKD)) |
189 |
|
|
DLF(JN) = REAL(2*JN-1,JPKD)/REAL(JN,JPKD) |
190 |
|
|
DLG(JN) = REAL(JN-1,JPKD)/REAL(JN,JPKD) |
191 |
|
|
DLH(JN) = SQRT(REAL(2*JN+1,JPKD)/REAL(2*JN,JPKD)) |
192 |
|
|
DLI(JN) = REAL(JN,JPKD) |
193 |
|
|
ENDDO |
194 |
|
|
|
195 |
|
|
!CALL GSTATS(1801,0) ! MPL 4.12.08 |
196 |
|
|
DO JGL=D%NLATLS(MYSETW),D%NLATLE(MYSETW) |
197 |
|
|
DLPOL(:,:) = 0.0_JPRB |
198 |
|
|
CALL SUPOL(R%NTMAX+1,DLRMU(JGL),DLPOL,DLA,DLB,DLC,DLD,DLE,DLF,DLG,DLH,DLI) |
199 |
|
|
INM = 0 |
200 |
|
|
IGLLOC = JGL - D%NLATLS(MYSETW) + 1 |
201 |
|
|
DO JM=0,R%NSMAX |
202 |
|
|
DO JN=R%NTMAX+1,JM,-1 |
203 |
|
|
INM = INM+1 |
204 |
|
|
ZPNMG(INM,IGLLOC) = REAL(DLPOL(JM,JN),JPKS) |
205 |
|
|
ENDDO |
206 |
|
|
ENDDO |
207 |
|
|
ENDDO |
208 |
|
|
!CALL GSTATS(1801,1) ! MPL 4.12.08 |
209 |
|
|
!CALL GSTATS(140,1) ! MPL 4.12.08 |
210 |
|
|
|
211 |
|
|
!CALL GSTATS(190,0) ! MPL 4.12.08 |
212 |
|
|
CALL SUTRLE(ZPNMG) |
213 |
|
|
!CALL GSTATS(190,1) ! MPL 4.12.08 |
214 |
|
|
|
215 |
|
|
ICOUNT = 0 |
216 |
|
|
DO JMLOC=1,D%NUMP |
217 |
|
|
IM = D%MYMS(JMLOC) |
218 |
|
|
DO JN=IM,R%NTMAX+2 |
219 |
|
|
ICOUNT = ICOUNT+1 |
220 |
|
|
ENDDO |
221 |
|
|
ENDDO |
222 |
|
|
|
223 |
|
|
ALLOCATE(F%REPSNM(ICOUNT)) |
224 |
|
|
IF (LLP2) WRITE(NOUT,9) 'F%REPSNM ',SIZE(F%REPSNM ),SHAPE(F%REPSNM ) |
225 |
|
|
|
226 |
|
|
ICOUNT = 0 |
227 |
|
|
DO JMLOC=1,D%NUMP |
228 |
|
|
IM = D%MYMS(JMLOC) |
229 |
|
|
DO JN=IM,R%NTMAX+2 |
230 |
|
|
ICOUNT = ICOUNT+1 |
231 |
|
|
F%REPSNM(ICOUNT) = REAL(SQRT(REAL(JN*JN-IM*IM,JPKD)/& |
232 |
|
|
&REAL(4*JN*JN-1,JPKD)),JPKS) |
233 |
|
|
ENDDO |
234 |
|
|
ENDDO |
235 |
|
|
|
236 |
|
|
ALLOCATE(F%RN(-1:R%NTMAX+3)) |
237 |
|
|
IF (LLP2) WRITE(NOUT,9) 'F%RN ',SIZE(F%RN ),SHAPE(F%RN ) |
238 |
|
|
ALLOCATE(F%RLAPIN(-1:R%NSMAX+2)) |
239 |
|
|
IF (LLP2) WRITE(NOUT,9) 'F%RLAPIN ',SIZE(F%RLAPIN ),SHAPE(F%RLAPIN ) |
240 |
|
|
ALLOCATE(F%NLTN(-1:R%NTMAX+3)) |
241 |
|
|
IF (LLP2) WRITE(NOUT,9) 'F%NLTN ',SIZE(F%NLTN ),SHAPE(F%NLTN ) |
242 |
|
|
|
243 |
|
|
DO JN=-1,R%NTMAX+3 |
244 |
|
|
F%RN(JN) = REAL(JN,JPRB) |
245 |
|
|
F%NLTN(JN) = R%NTMAX+2-JN |
246 |
|
|
ENDDO |
247 |
|
|
F%RLAPIN(:) = 0.0_JPRB |
248 |
|
|
F%RLAPIN(0) = 0._JPRB |
249 |
|
|
F%RLAPIN(-1) = 0.0_JPRB |
250 |
|
|
DO JN=1,R%NSMAX+2 |
251 |
|
|
F%RLAPIN(JN)=REAL(-(REAL(RA,JPKD)*REAL(RA,JPKD))/REAL(JN*(JN+1),JPKD),JPKS) |
252 |
|
|
ENDDO |
253 |
|
|
|
254 |
|
|
DEALLOCATE(ZPNMG) |
255 |
|
|
|
256 |
|
|
! ------------------------------------------------------------------ |
257 |
|
|
9 FORMAT(1X,'ARRAY ',A10,' ALLOCATED ',8I8) |
258 |
|
|
|
259 |
|
|
END SUBROUTINE SULEG |
260 |
|
|
END MODULE SULEG_MOD |