root/GCM/modelII/trunk/FFT36macDBL.f

Revision 10, 12.9 kB (checked in by mshopsin, 3 years ago)

model version 1.0.7

Line 
1c  *********************************************************************
2c  *********************************************************************
3c  **
4c  ** Model IImac
5c  ** Based on GCMII code for IBM RS/6000 computers created at GISS
6c  ** Modified to compile under Absoft Pro Fortran 6.2 for MacOS. 
7c  ** Based on MP008macC9, BA94C9 and MA94DC9
8c  **
9c  ** CHANGE HISTORY:
10c  **
11c  ** 07/27/99 First Successful Compile! (DCH)
12c  ** 12/21/00 converted for new model (MFS)
13c  **
14c  ** NOTES:
15c  **
16c  *********************************************************************
17c  *********************************************************************
18
19      SUBROUTINE FRTR0 (KM)
20c  ** INITIALIZATION ENTRY TO CALCULATE SIN VALUES AND CHECK THAT KM=36
21      IMPLICIT REAL*8 (A-H,O-Z)
22      COMMON/FCOM/BYKM,BYKMH,BYKM2,SIN10,SIN20,SIN30,SIN40,SIN50,SIN60,
23     *  SIN70,SIN80
24c     REAL*8 TWOPI/6.283185307179586477/
25      DATA   TWOPI/6.283185307179586477/
26      IF(KM.NE.36) GO TO 220
27      BYKM=1./KM
28      BYKMH=2./KM
29      BYKM2=1./(2.*KM)
30      SIN10=DSIN(TWOPI/36.)
31      SIN20=DSIN(TWOPI/18.)
32      SIN30=1.D0/2.
33      SIN40=DSIN(TWOPI/9.)
34      SIN50=DCOS(TWOPI/9.)
35      SIN60=DSQRT(3.D0)/2.
36      SIN70=DCOS(TWOPI/18.)
37      SIN80=DCOS(TWOPI/36.)
38c  BL SIN10=DSIN(TWOPI/36.)
39c  BL SIN20=DSIN(TWOPI/18.)
40c  BL SIN30=1.D0/2.
41c  BL SIN40=DSIN(TWOPI/9.)
42c  BL SIN50=DCOS(TWOPI/9.)
43c  BL SIN60=DSQRT(3.D0)/2.
44c  BL SIN70=DCOS(TWOPI/18.)
45c  BL SIN80=DCOS(TWOPI/36.)
46c     SIN10=SIN(TWOPI/36.)
47c     SIN20=SIN(TWOPI/18.)
48c     SIN30=.5
49c     SIN40=SIN(TWOPI/9.)
50c     SIN50=COS(TWOPI/9.)
51c     SIN60=SQRT(3.)/2.
52c     SIN70=COS(TWOPI/18.)
53c     SIN80=COS(TWOPI/36.)
54      RETURN
55  220 WRITE (6,901) KM
56      STOP
57  901 FORMAT ('0THIS FOURT SUBROUTINE NOT SUITED FOR KM = ',I8)
58      END
59      SUBROUTINE FRTR (F)
60c  ** THIS SUBROUTINE PERFORMS A FOURIER ANALYSIS ON THE ONE DIMENSIONAL
61c  ** ARRAY F WHICH MUST BE DIMENSIONED 36.  IT RETURNS IN F THE ENERGY
62c  ** ASSOCIATED WITH EACH WAVE NUMBER.  UPON ENTERING THIS ROUTINE,
63c  ** THE TOTAL ENERGY IS
64c  **   .5*SUM(F(K)*F(K))
65c  ** WITH THE SUM BEING TAKEN OVER ALL K FROM 1 TO 36.  UPON LEAVING
66c  ** THIS ROUTINE, THE TOTAL ENERGY IS
67c  **   SUM(F(N+1))
68c  ** WITH THE SUM BEING TAKEN OVER ALL WAVE NUMBERS FROM 0 TO 18.
69      IMPLICIT REAL*8 (A-H,O-Z)
70      COMMON/FCOM/BYKM,BYKMH,BYKM2,SIN10,SIN20,SIN30,SIN40,SIN50,SIN60,
71     *  SIN70,SIN80
72      DIMENSION F(36)
73   10 CC00=F(12)+F(24)+F(36)
74      CC01=F(36)-(F(12)+F(24))*SIN30
75      CC10=F(1)+F(13)+F(25)
76      CC11=F(1)*SIN80-F(13)*SIN40-F(25)*SIN20
77      CC20=F(2)+F(14)+F(26)
78      CC21=F(2)*SIN70-F(14)*SIN50-F(26)*SIN10
79      CC30=F(3)+F(15)+F(27)
80      CC31=(F(3)-F(15))*SIN60
81      CC40=F(4)+F(16)+F(28)
82      CC41=F(4)*SIN50-F(16)*SIN70+F(28)*SIN10
83      CC50=F(5)+F(17)+F(29)
84      CC51=F(5)*SIN40-F(17)*SIN80+F(29)*SIN20
85      CC60=F(6)+F(18)+F(30)
86      CC61=(F(6)+F(30))*SIN30-F(18)
87      CC70=F(7)+F(19)+F(31)
88      CC71=F(7)*SIN20-F(19)*SIN80+F(31)*SIN40
89      CC80=F(8)+F(20)+F(32)
90      CC81=F(8)*SIN10-F(20)*SIN70+F(32)*SIN50
91      CC90=F(9)+F(21)+F(33)
92      CC91=(F(33)-F(21))*SIN60
93      CCA0=F(10)+F(22)+F(34)
94      CCA1=F(34)*SIN70-F(10)*SIN10-F(22)*SIN50
95      CCB0=F(11)+F(23)+F(35)
96      CCB1=F(35)*SIN80-F(11)*SIN20-F(23)*SIN40
97      SC01=(F(12)-F(24))*SIN60
98      SC11=F(1)*SIN10+F(13)*SIN50-F(25)*SIN70
99      SC21=F(2)*SIN20+F(14)*SIN40-F(26)*SIN80
100      SC31=(F(3)+F(15))*SIN30-F(27)
101      SC41=F(4)*SIN40+F(16)*SIN20-F(28)*SIN80
102      SC51=F(5)*SIN50+F(17)*SIN10-F(29)*SIN70
103      SC61=(F(6)-F(30))*SIN60
104      SC71=F(7)*SIN70-F(19)*SIN10-F(31)*SIN50
105      SC81=F(8)*SIN80-F(20)*SIN20-F(32)*SIN40
106      SC91=F(9)-(F(21)+F(33))*SIN30
107      SCA1=F(10)*SIN80-F(22)*SIN40-F(34)*SIN20
108      SCB1=F(11)*SIN70-F(23)*SIN50-F(35)*SIN10
109c  ** CALCULATE EXPRESSIONS SUMMED BY INCREMENTS OF 4
110      C400=CC00+CC40+CC80
111      C401=CC01+CC41+CC81
112      C403=CC00-(CC40+CC80)*SIN30
113      C402=(CC01-(CC41+CC81)*SIN30)+((SC41-SC81)*SIN60)
114      C404=(CC01-(CC41+CC81)*SIN30)-((SC41-SC81)*SIN60)
115      C410=CC10+CC50+CC90
116      C411=CC11+CC51+CC91
117      C413=(CC10-CC50)*SIN60
118      C412=((CC11-CC51)*SIN60)+((SC11+SC51)*SIN30-SC91)
119      C414=((CC11-CC51)*SIN60)-((SC11+SC51)*SIN30-SC91)
120      C420=CC20+CC60+CCA0
121      C421=CC21+CC61+CCA1
122      C423=(CC20+CCA0)*SIN30-CC60
123      C422=((CC21+CCA1)*SIN30-CC61)+((SC21-SCA1)*SIN60)
124      C424=((CC21+CCA1)*SIN30-CC61)-((SC21-SCA1)*SIN60)
125      C430=CC30+CC70+CCB0
126      C431=CC31+CC71+CCB1
127      C433=(CCB0-CC70)*SIN60
128      C432=((CCB1-CC71)*SIN60)+(SC31-(SC71+SCB1)*SIN30)
129      C434=((CCB1-CC71)*SIN60)-(SC31-(SC71+SCB1)*SIN30)
130      S401=SC01+SC41+SC81
131      S403=(CC40-CC80)*SIN60
132      S402=((CC41-CC81)*SIN60)+((SC41+SC81)*SIN30-SC01)
133      S404=((CC41-CC81)*SIN60)-((SC41+SC81)*SIN30-SC01)
134      S411=SC11+SC51+SC91
135      S413=(CC10+CC50)*SIN30-CC90
136      S412=((CC11+CC51)*SIN30-CC91)+((SC51-SC11)*SIN60)
137      S414=((CC11+CC51)*SIN30-CC91)-((SC51-SC11)*SIN60)
138      S421=SC21+SC61+SCA1
139      S423=(CC20-CCA0)*SIN60
140      S422=((CC21-CCA1)*SIN60)+(SC61-(SC21+SCA1)*SIN30)
141      S424=((CC21-CCA1)*SIN60)-(SC61-(SC21+SCA1)*SIN30)
142      S431=SC31+SC71+SCB1
143      S433=CC30-(CC70+CCB0)*SIN30
144      S432=(CC31-(CC71+CCB1)*SIN30)+((SC71-SCB1)*SIN60)
145      S434=(CC31-(CC71+CCB1)*SIN30)-((SC71-SCB1)*SIN60)
146c  ** CALCULATE EXPRESSIONS SUMMED BY INCREMENTS OF 2
147      C200=C400+C420
148      C201=C401+C421
149      C202=C402+C422
150      C203=C403+C423
151      C204=C404+C424
152      C205=C404-C424
153      C206=C403-C423
154      C207=C402-C422
155      C208=C401-C421
156c     C209=C400-C420
157      C210=C410+C430
158      C211=C411+C431
159      C212=C412+C432
160      C213=C413+C433
161      C214=C414+C434
162      C215=S414-S434
163      C216=S413-S433
164      C217=S412-S432
165      C218=S411-S431
166c     C219=0
167c     S200=0
168      S201=S401+S421
169      S202=S402+S422
170      S203=S403+S423
171      S204=S404+S424
172      S205=S424-S404
173      S206=S423-S403
174      S207=S422-S402
175      S208=S421-S401
176c     S209=0
177c     S210=0
178      S211=S411+S431
179      S212=S412+S432
180      S213=S413+S433
181      S214=S414+S434
182      S215=C414-C434
183      S216=C413-C433
184      S217=C412-C432
185      S218=C411-C431
186c     S219=C410-C430
187c  ** CALCULATE THE SQUARE OF THE MAGNITUDE OF G(1,N)+I*G(2,N)
188   20 F(1)=(C200+C210)*(C200+C210)*BYKM2
189      F(2)=((C201+C211)*(C201+C211)+(S201+S211)*(S201+S211))*BYKM
190      F(3)=((C202+C212)*(C202+C212)+(S202+S212)*(S202+S212))*BYKM
191      F(4)=((C203+C213)*(C203+C213)+(S203+S213)*(S203+S213))*BYKM
192      F(5)=((C204+C214)*(C204+C214)+(S204+S214)*(S204+S214))*BYKM
193      F(6)=((C205+C215)*(C205+C215)+(S205+S215)*(S205+S215))*BYKM
194      F(7)=((C206+C216)*(C206+C216)+(S206+S216)*(S206+S216))*BYKM
195      F(8)=((C207+C217)*(C207+C217)+(S207+S217)*(S207+S217))*BYKM
196      F(9)=((C208+C218)*(C208+C218)+(S208+S218)*(S208+S218))*BYKM
197      F(10)=((C400-C420)*(C400-C420)+(C410-C430)*(C410-C430))*BYKM
198      F(11)=((C208-C218)*(C208-C218)+(S218-S208)*(S218-S208))*BYKM
199      F(12)=((C207-C217)*(C207-C217)+(S217-S207)*(S217-S207))*BYKM
200      F(13)=((C206-C216)*(C206-C216)+(S216-S206)*(S216-S206))*BYKM
201      F(14)=((C205-C215)*(C205-C215)+(S215-S205)*(S215-S205))*BYKM
202      F(15)=((C204-C214)*(C204-C214)+(S214-S204)*(S214-S204))*BYKM
203      F(16)=((C203-C213)*(C203-C213)+(S213-S203)*(S213-S203))*BYKM
204      F(17)=((C202-C212)*(C202-C212)+(S212-S202)*(S212-S202))*BYKM
205      F(18)=((C201-C211)*(C201-C211)+(S211-S201)*(S211-S201))*BYKM
206      F(19)=(C200-C210)*(C200-C210)*BYKM2
207      RETURN
208      END
209      SUBROUTINE GETAN (F,G)
210c  ** GETAN RETRIEVES THE FOURIER COEFFICIENTS CONTAINED IN AN
211c  ** ARRAY G DIMENSIONED 2 BY 19 AND DEFINED BY
212c  **   G(1,N+1)+I*G(2,N+1)=SUM(F(K)*EXP(-2*PI*I*N*K/KM))/KMH
213c  ** WITH THE SUM TAKEN OVER ALL K FROM 1 TO KM.  KMH = KM FOR N = 0
214c  ** OR 18, OTHERWISE KMH = KM/2.  THE INTERNAL NOTATION CPQN MEANS
215c  **   CPQN = SUM(F(K)*COS(2*PI*N*K/KM))
216c  ** WITH THE SUM BEING TAKEN OVER ALL K FROM 1 TO KM WHICH ARE EQUAL
217c  ** TO Q MODULO(P).  SPQN IS THE SAME BUT WITH COS REPLACED BY SIN.
218c  ** THE NOTATION A=10, B=11, ETC. IS USED FOR  P, Q AND N.
219      IMPLICIT REAL*8 (A-H,O-Z)
220      COMMON/FCOM/BYKM,BYKMH,BYKM2,SIN10,SIN20,SIN30,SIN40,SIN50,SIN60,
221     *  SIN70,SIN80
222      DIMENSION F(36),G(2,19)
223c  ** CALCULATE EXPRESSIONS SUMMED BY INCREMENTS OF 12
224   10 CC00=F(12)+F(24)+F(36)
225      CC01=F(36)-(F(12)+F(24))*SIN30
226      CC10=F(1)+F(13)+F(25)
227      CC11=F(1)*SIN80-F(13)*SIN40-F(25)*SIN20
228      CC20=F(2)+F(14)+F(26)
229      CC21=F(2)*SIN70-F(14)*SIN50-F(26)*SIN10
230      CC30=F(3)+F(15)+F(27)
231      CC31=(F(3)-F(15))*SIN60
232      CC40=F(4)+F(16)+F(28)
233      CC41=F(4)*SIN50-F(16)*SIN70+F(28)*SIN10
234      CC50=F(5)+F(17)+F(29)
235      CC51=F(5)*SIN40-F(17)*SIN80+F(29)*SIN20
236      CC60=F(6)+F(18)+F(30)
237      CC61=(F(6)+F(30))*SIN30-F(18)
238      CC70=F(7)+F(19)+F(31)
239      CC71=F(7)*SIN20-F(19)*SIN80+F(31)*SIN40
240      CC80=F(8)+F(20)+F(32)
241      CC81=F(8)*SIN10-F(20)*SIN70+F(32)*SIN50
242      CC90=F(9)+F(21)+F(33)
243      CC91=(F(33)-F(21))*SIN60
244      CCA0=F(10)+F(22)+F(34)
245      CCA1=F(34)*SIN70-F(10)*SIN10-F(22)*SIN50
246      CCB0=F(11)+F(23)+F(35)
247      CCB1=F(35)*SIN80-F(11)*SIN20-F(23)*SIN40
248      SC01=(F(12)-F(24))*SIN60
249      SC11=F(1)*SIN10+F(13)*SIN50-F(25)*SIN70
250      SC21=F(2)*SIN20+F(14)*SIN40-F(26)*SIN80
251      SC31=(F(3)+F(15))*SIN30-F(27)
252      SC41=F(4)*SIN40+F(16)*SIN20-F(28)*SIN80
253      SC51=F(5)*SIN50+F(17)*SIN10-F(29)*SIN70
254      SC61=(F(6)-F(30))*SIN60
255      SC71=F(7)*SIN70-F(19)*SIN10-F(31)*SIN50
256      SC81=F(8)*SIN80-F(20)*SIN20-F(32)*SIN40
257      SC91=F(9)-(F(21)+F(33))*SIN30
258      SCA1=F(10)*SIN80-F(22)*SIN40-F(34)*SIN20
259      SCB1=F(11)*SIN70-F(23)*SIN50-F(35)*SIN10
260c  ** CALCULATE EXPRESSIONS SUMMED BY INCREMENTS OF 4
261      C400=CC00+CC40+CC80
262      C401=CC01+CC41+CC81
263      C403=CC00-(CC40+CC80)*SIN30
264      C402=(CC01-(CC41+CC81)*SIN30)+((SC41-SC81)*SIN60)
265      C404=(CC01-(CC41+CC81)*SIN30)-((SC41-SC81)*SIN60)
266      C410=CC10+CC50+CC90
267      C411=CC11+CC51+CC91
268      C413=(CC10-CC50)*SIN60
269      C412=((CC11-CC51)*SIN60)+((SC11+SC51)*SIN30-SC91)
270      C414=((CC11-CC51)*SIN60)-((SC11+SC51)*SIN30-SC91)
271      C420=CC20+CC60+CCA0
272      C421=CC21+CC61+CCA1
273      C423=(CC20+CCA0)*SIN30-CC60
274      C422=((CC21+CCA1)*SIN30-CC61)+((SC21-SCA1)*SIN60)
275      C424=((CC21+CCA1)*SIN30-CC61)-((SC21-SCA1)*SIN60)
276      C430=CC30+CC70+CCB0
277      C431=CC31+CC71+CCB1
278      C433=(CCB0-CC70)*SIN60
279      C432=((CCB1-CC71)*SIN60)+(SC31-(SC71+SCB1)*SIN30)
280      C434=((CCB1-CC71)*SIN60)-(SC31-(SC71+SCB1)*SIN30)
281      S401=SC01+SC41+SC81
282      S403=(CC40-CC80)*SIN60
283      S402=((CC41-CC81)*SIN60)+((SC41+SC81)*SIN30-SC01)
284      S404=((CC41-CC81)*SIN60)-((SC41+SC81)*SIN30-SC01)
285      S411=SC11+SC51+SC91
286      S413=(CC10+CC50)*SIN30-CC90
287      S412=((CC11+CC51)*SIN30-CC91)+((SC51-SC11)*SIN60)
288      S414=((CC11+CC51)*SIN30-CC91)-((SC51-SC11)*SIN60)
289      S421=SC21+SC61+SCA1
290      S423=(CC20-CCA0)*SIN60
291      S422=((CC21-CCA1)*SIN60)+(SC61-(SC21+SCA1)*SIN30)
292      S424=((CC21-CCA1)*SIN60)-(SC61-(SC21+SCA1)*SIN30)
293      S431=SC31+SC71+SCB1
294      S433=CC30-(CC70+CCB0)*SIN30
295      S432=(CC31-(CC71+CCB1)*SIN30)+((SC71-SCB1)*SIN60)
296      S434=(CC31-(CC71+CCB1)*SIN30)-((SC71-SCB1)*SIN60)
297c  ** CALCULATE EXPRESSIONS SUMMED BY INCREMENTS OF 2
298      C200=C400+C420
299      C201=C401+C421
300      C202=C402+C422
301      C203=C403+C423
302      C204=C404+C424
303      C205=C404-C424
304      C206=C403-C423
305      C207=C402-C422
306      C208=C401-C421
307c     C209=C400-C420
308      C210=C410+C430
309      C211=C411+C431
310      C212=C412+C432
311      C213=C413+C433
312      C214=C414+C434
313      C215=S414-S434
314      C216=S413-S433
315      C217=S412-S432
316      C218=S411-S431
317c     C219=0
318c     S200=0
319      S201=S401+S421
320      S202=S402+S422
321      S203=S403+S423
322      S204=S404+S424
323      S205=S424-S404
324      S206=S423-S403
325      S207=S422-S402
326      S208=S421-S401
327c     S209=0
328c     S210=0
329      S211=S411+S431
330      S212=S412+S432
331      S213=S413+S433
332      S214=S414+S434
333      S215=C414-C434
334      S216=C413-C433
335      S217=C412-C432
336      S218=C411-C431
337c     S219=C410-C430
338c  ** CALCULATE FINAL COEFFICIENTS OF FOURIER EXPANSION
339      G(1,1)=(C200+C210)*BYKM
340      G(1,2)=(C201+C211)*BYKMH
341      G(1,3)=(C202+C212)*BYKMH
342      G(1,4)=(C203+C213)*BYKMH
343      G(1,5)=(C204+C214)*BYKMH
344      G(1,6)=(C205+C215)*BYKMH
345      G(1,7)=(C206+C216)*BYKMH
346      G(1,8)=(C207+C217)*BYKMH
347      G(1,9)=(C208+C218)*BYKMH
348      G(1,10)=(C400-C420)*BYKMH
349      G(1,11)=(C208-C218)*BYKMH
350      G(1,12)=(C207-C217)*BYKMH
351      G(1,13)=(C206-C216)*BYKMH
352      G(1,14)=(C205-C215)*BYKMH
353      G(1,15)=(C204-C214)*BYKMH
354      G(1,16)=(C203-C213)*BYKMH
355      G(1,17)=(C202-C212)*BYKMH
356      G(1,18)=(C201-C211)*BYKMH
357      G(1,19)=(C200-C210)*BYKM
358      G(2,1)=0.
359      G(2,2)=(S201+S211)*BYKMH
360      G(2,3)=(S202+S212)*BYKMH
361      G(2,4)=(S203+S213)*BYKMH
362      G(2,5)=(S204+S214)*BYKMH
363      G(2,6)=(S205+S215)*BYKMH
364      G(2,7)=(S206+S216)*BYKMH
365      G(2,8)=(S207+S217)*BYKMH
366      G(2,9)=(S208+S218)*BYKMH
367      G(2,10)=(C410-C430)*BYKMH
368      G(2,11)=(S218-S208)*BYKMH
369      G(2,12)=(S217-S207)*BYKMH
370      G(2,13)=(S216-S206)*BYKMH
371      G(2,14)=(S215-S205)*BYKMH
372      G(2,15)=(S214-S204)*BYKMH
373      G(2,16)=(S213-S203)*BYKMH
374      G(2,17)=(S212-S202)*BYKMH
375      G(2,18)=(S211-S201)*BYKMH
376      G(2,19)=0.
377      RETURN
378      END
Note: See TracBrowser for help on using the browser.