Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

swr_toon.F @ 317

Last change on this file since 317 was 38, checked in by emillour, 14 years ago
Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs. EM
File size: 17.5 KB

Line
1	SUBROUTINE SWR_TOON ( KDLON, KFLEV, KNU
2	S , aerosol,QVISsQREF3d,omegaVIS3d,gVIS3d
3	& , albedo,PDSIG,PPSOL,PRMU,PSEC
4	S , PFD,PFU )
5
6	IMPLICIT NONE
7	C
8	#include "dimensions.h"
9	#include "dimphys.h"
10	#include "dimradmars.h"
11	#include "callkeys.h"
12
13	#include "yomaer.h"
14	#include "yomlw.h"
15
16	C
17	C SWR - Continuum scattering computations
18	C
19	C PURPOSE.
20	C --------
21	C Computes the reflectivity and transmissivity in case oF
22	C Continuum scattering
23	c F. Forget (1999)
24	c
25	c Modified by Tran The Trung, using radiative transfer code
26	c of Toon 1981.
27	C
28	C IMPLICIT ARGUMENTS :
29	C --------------------
30	C
31	C ==== INPUTS ===
32	c
33	c KDLON : number of horizontal grid points
34	c KFLEV : number of vertical layers
35	c KNU : Solar band # (1 or 2)
36	c aerosol aerosol extinction optical depth
37	c at reference wavelength "longrefvis" set
38	c in dimradmars.h , in each layer, for one of
39	c the "naerkind" kind of aerosol optical properties.
40	c albedo hemispheric surface albedo
41	c albedo (i,1) : mean albedo for solar band#1
42	c (see below)
43	c albedo (i,2) : mean albedo for solar band#2
44	c (see below)
45	c PDSIG layer thickness in sigma coordinates
46	c PPSOL Surface pressure (Pa)
47	c PRMU: cos of solar zenith angle (=1 when sun at zenith)
48	c (CORRECTED for high zenith angle (atmosphere), unlike mu0)
49	c PSEC =1./PRMU
50
51	C ==== OUTPUTS ===
52	c
53	c PFD : downward flux in spectral band #INU in a given mesh
54	c (normalized to the total incident flux at the top of the atmosphere)
55	c PFU : upward flux in specatral band #INU in a given mesh
56	c (normalized to the total incident flux at the top of the atmosphere)
57	C
58	C
59	C METHOD.
60	C -------
61	C
62	C Computes continuum fluxes corresponding to aerosoL
63	C Or/and rayleigh scattering (no molecular gas absorption)
64	C
65	C-----------------------------------------------------------------------
66	C
67	C
68	C-----------------------------------------------------------------------
69	C
70
71	C ARGUMENTS
72	C ---------
73	INTEGER KDLON, KFLEV, KNU
74	REAL aerosol(NDLO2,KFLEV,naerkind), albedo(NDLO2,2),
75	S PDSIG(NDLO2,KFLEV),PSEC(NDLO2)
76
77	REAL QVISsQREF3d(NDLO2,KFLEV,nsun,naerkind)
78	REAL omegaVIS3d(NDLO2,KFLEV,nsun,naerkind)
79	REAL gVIS3d(NDLO2,KFLEV,nsun,naerkind)
80
81	REAL PPSOL(NDLO2)
82	REAL PFD(NDLO2,KFLEV+1),PFU(NDLO2,KFLEV+1)
83	REAL PRMU(NDLO2)
84
85	C LOCAL ARRAYS
86	C ------------
87
88	INTEGER jk,jl,jae
89	REAL ZTRAY, ZRATIO,ZGAR, ZFF
90	REAL ZCGAZ(NDLO2,NFLEV),ZPIZAZ(NDLO2,NFLEV),ZTAUAZ(NDLO2,NFLEV)
91	REAL ZRAYL(NDLON)
92
93	c Part added by Tran The Trung
94	c inputs to gfluxv
95	REAL*8 DTDEL(nlaylte), WDEL(nlaylte), CDEL(nlaylte)
96	c outputs of gfluxv
97	REAL*8 FP(nlaylte+1), FM(nlaylte+1)
98	c normalization of top downward flux
99	REAL*8 norm
100	c End part added by Tran The Trung
101
102	c Function
103	c --------
104	real CVMGT
105
106	c Computing TOTAL single scattering parameters by adding
107	c properties of all the NAERKIND kind of aerosols
108
109	DO JK = 1 , nlaylte
110	DO JL = 1 , KDLON
111	ZCGAZ(JL,JK) = 0.
112	ZPIZAZ(JL,JK) = 0.
113	ZTAUAZ(JL,JK) = 0.
114	END DO
115	DO 106 JAE=1,naerkind
116	DO 105 JL = 1 , KDLON
117	c Mean Extinction optical depth in the spectral band
118	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
119	ZTAUAZ(JL,JK)=ZTAUAZ(JL,JK)
120	S +aerosol(JL,JK,JAE)*QVISsQREF3d(JL,JK,KNU,JAE)
121	c Single scattering albedo
122	c ~~~~~~~~~~~~~~~~~~~~~~~~
123	ZPIZAZ(JL,JK)=ZPIZAZ(JL,JK)+aerosol(JL,JK,JAE)*
124	S QVISsQREF3d(JL,JK,KNU,JAE)*
125	& omegaVIS3d(JL,JK,KNU,JAE)
126	c Assymetry factor
127	c ~~~~~~~~~~~~~~~~
128	ZCGAZ(JL,JK) = ZCGAZ(JL,JK) +aerosol(JL,JK,JAE)*
129	S QVISsQREF3d(JL,JK,KNU,JAE)*
130	& omegaVIS3d(JL,JK,KNU,JAE)*gVIS3d(JL,JK,KNU,JAE)
131	105 CONTINUE
132	106 CONTINUE
133	END DO
134	C
135	DO JK = 1 , nlaylte
136	DO JL = 1 , KDLON
137	ZCGAZ(JL,JK) = CVMGT( 0., ZCGAZ(JL,JK) / ZPIZAZ(JL,JK),
138	S (ZPIZAZ(JL,JK).EQ.0) )
139	ZPIZAZ(JL,JK) = CVMGT( 1., ZPIZAZ(JL,JK) / ZTAUAZ(JL,JK),
140	S (ZTAUAZ(JL,JK).EQ.0) )
141	END DO
142	END DO
143
144	C --------------------------------
145	C INCLUDING RAYLEIGH SCATERRING
146	C -------------------------------
147	if (rayleigh) then
148
149	call swrayleigh(kdlon,knu,ppsol,prmu,ZRAYL)
150
151	c Modifying mean aerosol parameters to account rayleigh scat by gas:
152
153	DO JK = 1 , nlaylte
154	DO JL = 1 , KDLON
155	c Rayleigh opacity in each layer :
156	ZTRAY = ZRAYL(JL) * PDSIG(JL,JK)
157	c ratio Tau(rayleigh) / Tau (total)
158	ZRATIO = ZTRAY / (ZTRAY + ZTAUAZ(JL,JK))
159	ZGAR = ZCGAZ(JL,JK)
160	ZFF = ZGAR * ZGAR
161	ZTAUAZ(JL,JK)=ZTRAY+ZTAUAZ(JL,JK)(1.-ZPIZAZ(JL,JK)ZFF)
162	ZCGAZ(JL,JK) = ZGAR * (1. - ZRATIO) / (1. + ZGAR)
163	ZPIZAZ(JL,JK) =ZRATIO+(1.-ZRATIO)ZPIZAZ(JL,JK)(1.-ZFF)
164	S / (1. -ZPIZAZ(JL,JK) * ZFF)
165	END DO
166	END DO
167	end if
168
169	c Part added by Tran The Trung
170	c new radiative transfer
171
172	do JL = 1, KDLON
173
174	c assign temporary inputs
175	do JK = 1, nlaylte
176	jae = nlaylte+1-JK
177	DTDEL(JK) = real(ZTAUAZ(JL,jae),8)
178	WDEL(JK) = real(ZPIZAZ(JL,jae),8)
179	CDEL(JK) = real(ZCGAZ(JL,jae),8)
180	end do
181
182	c call gfluxv
183	call gfluxv(nlaylte, DTDEL,WDEL,CDEL,
184	S real(PRMU(JL),8),
185	S real(albedo(JL,KNU),8),
186	S FP,FM)
187
188	c assign output
189	norm = FM(1)
190	c here we can have a check of norm not being 0.0
191	c however it would never happen in practice,
192	c so we can comment out
193	c if (norm .gt. 0.0) then
194	do JK = 1, nlaylte+1
195	jae = nlaylte+2-JK
196	PFU(JL,JK) = sunfr(KNU)*real(FP(jae)/norm,4)
197	PFD(JL,JK) = sunfr(KNU)*real(FM(jae)/norm,4)
198	end do
199	c elseif (norm .eq. 0.0) then
200	c do JK = 1, nlaylte+1
201	c PFU(JL,JK) = 0.0
202	c PFD(JL,JK) = 0.0
203	c end do
204	c else
205	c stop "Error: top downward visible flux is negative!"
206	c end if
207
208	end do
209	c End part added by Tran The Trung
210
211	RETURN
212	END
213
214	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
215
216	SUBROUTINE GFLUXV(NAYER,DTDEL,WDEL,CDEL,UBAR0,RSF
217	& ,FP,FM)
218	IMPLICIT NONE
219	C THIS SUBROUTINE TAKES THE OPTICAL CONSTANTS AND BOUNDARY CONDITONS
220	C FOR THE VISIBLE FLUX AT ONE WAVELENGTH AND SOLVES FOR THE FLUXES AT
221	C THE LEVELS. THIS VERSION IS SET UP TO WORK WITH LAYER OPTICAL DEPTHS
222	C MEASURED FROM THE TOP OF EACH LAYER. (DTAU) TOP OF EACH LAYER HAS
223	C OPTICAL DEPTH TAU(N). IN THIS SUB LEVEL N IS ABOVE LAYER N. THAT IS LAYER N
224	C HAS LEVEL N ON TOP AND LEVEL N+1 ON BOTTOM. OPTICAL DEPTH INCREASES
225	C FROM TOP TO BOTTOM. SEE C.P. MCKAY, TGM NOTES.
226	C THIS SUBROUTINE DIFFERS FROM ITS IR CONTERPART IN THAT HERE WE SOLVE FOR
227	C THE FLUXES DIRECTLY USING THE GENERALIZED NOTATION OF MEADOR AND WEAVOR
228	C J.A.S., 37, 630-642, 1980.
229	C THE TRI-DIAGONAL MATRIX SOLVER IS DSOLVER AND IS DOUBLE PRECISION SO MANY
230	C VARIABLES ARE PASSED AS SINGLE THEN BECOME DOUBLE IN DSOLVER
231
232	C THIS VERSION HAS BEEN MODIFIED BY TRAN THE TRUNG WITH:
233	C 1. Simplified input & output for swr.F subroutine in LMDZ.MARS gcm model
234	C 2. Use delta function to modify optical properties
235	C 3. Use delta-eddington G1, G2, G3 parameters
236	C 4. Optimized for speed
237
238	C INPUTS:
239	INTEGER NAYER
240	c NAYER = number of layer
241	c first layer is at top
242	c last layer is at bottom
243	REAL*8 DTDEL(NAYER), WDEL(NAYER), CDEL(NAYER)
244	c DTDEL = optical thickness of layer
245	c WDEL = single scattering of layer
246	c CDEL = assymetry parameter
247	REAL*8 UBAR0, RSF
248	c UBAR0 = absolute value of cosine of solar zenith angle
249	c RSF = surface albedo
250	C OUTPUTS:
251	REAL*8 FP(NAYER+1), FM(NAYER+1)
252	c FP = flux up
253	c FM = flux down
254	C PRIVATES:
255	INTEGER J,NL,NLEV
256	!!!! AS+JBM 03/2010 BUG BUG si trop niveaux verticaux (LES)
257	!!!! ET PAS BESOIN DE HARDWIRE SALE ICI !
258	!!!! CORRIGER CE BUG AMELIORE EFFICACITE ET FLEXIBILITE
259	!! PARAMETER (NL=201)
260	!! C THIS VALUE (201) MUST BE .GE. 2*NAYER
261	REAL*8 BSURF,AP,AM,DENOM,EM,EP,G4
262	!! REAL*8 W0(NL), COSBAR(NL), DTAU(NL), TAU(NL)
263	!! REAL*8 LAMDA(NL),XK1(NL),XK2(NL)
264	!! REAL*8 G1(NL),G2(NL),G3(NL)
265	!! REAL*8 GAMA(NL),CP(NL),CM(NL),CPM1(NL),CMM1(NL)
266	!! REAL*8 E1(NL),E2(NL),E3(NL),E4(NL)
267	REAL8 W0(2NAYER), COSBAR(2NAYER), DTAU(2NAYER), TAU(2*NAYER)
268	REAL8 LAMDA(2NAYER),XK1(2NAYER),XK2(2NAYER)
269	REAL8 G1(2NAYER),G2(2NAYER),G3(2NAYER)
270	REAL8 GAMA(2NAYER),CP(2NAYER),CM(2NAYER),CPM1(2*NAYER)
271	REAL8 CMM1(2NAYER)
272	REAL8 E1(2NAYER),E2(2NAYER),E3(2NAYER),E4(2*NAYER)
273
274	NL = 2*NAYER !!! AS+JBM 03/2010
275	NLEV = NAYER+1
276
277	C TURN ON THE DELTA-FUNCTION IF REQUIRED HERE
278	c TAU(1) = 0.0
279	c DO J=1,NAYER
280	c W0(J)=WDEL(J)
281	c COSBAR(J)=CDEL(J)
282	c DTAU(J)=DTDEL(J)
283	c TAU(J+1)=TAU(J)+DTAU(J)
284	c END DO
285	C FOR THE DELTA FUNCTION HERE...
286	TAU(1) = 0.0
287	DO J=1,NAYER
288	COSBAR(J)=CDEL(J)/(1.+CDEL(J))
289	W0(J)=1.-WDEL(J)CDEL(J)*2
290	DTAU(J)=DTDEL(J)*W0(J)
291	W0(J)=WDEL(J)(1.-CDEL(J)*2)/W0(J)
292	TAU(J+1)=TAU(J)+DTAU(J)
293	END DO
294
295	c Optimization, this is the major speed gain
296	TAU(1) = 1.0
297	do J = 2, NAYER+1
298	TAU(J) = EXP(-TAU(J)/UBAR0)
299	end do
300	BSURF = RSFUBAR0TAU(NLEV)
301	C WE GO WITH THE HEMISPHERIC CONSTANT APPROACH
302	C AS DEFINED BY M&W - THIS IS THE WAY THE IR IS DONE
303	DO J=1,NAYER
304	c Optimization: ALPHA not used
305	c ALPHA(J)=SQRT( (1.-W0(J))/(1.-W0(J)*COSBAR(J)) )
306	C SET OF CONSTANTS DETERMINED BY DOM
307	c G1(J)= (SQRT(3.)0.5)(2. - W0(J)*(1.+COSBAR(J)))
308	c G2(J)= (SQRT(3.)W0(J)0.5)*(1.-COSBAR(J))
309	c G3(J)=0.5(1.-SQRT(3.)COSBAR(J)*UBAR0)
310	c We use delta-Eddington instead
311	G1(J)=0.25(7.-W0(j)(4.+3*cosbar(j)))
312	G2(J)=-0.25(1.-W0(j)(4.-3*cosbar(j)))
313	G3(J)=0.5(1.-SQRT(3.)COSBAR(J)*UBAR0)
314	LAMDA(J)=SQRT(G1(J)2 - G2(J)2)
315	GAMA(J)=(G1(J)-LAMDA(J))/G2(J)
316	END DO
317
318	DO J=1,NAYER
319	G4=1.-G3(J)
320	DENOM=LAMDA(J)2 - 1./UBAR02
321	C NOTE THAT THE ALGORITHM DONOT ACCEPT UBAR0 .eq. 0
322	C THERE IS A POTENTIAL PROBLEM HERE IF W0=0 AND UBARV=UBAR0
323	C THEN DENOM WILL VANISH. THIS ONLY HAPPENS PHYSICALLY WHEN
324	C THE SCATTERING GOES TO ZERO
325	C PREVENT THIS WITH AN IF STATEMENT
326	IF ( DENOM .EQ. 0.) THEN
327	DENOM=1.E-6
328	END IF
329	DENOM = W0(J)/DENOM
330	AM=DENOM(G4 (G1(J)+1./UBAR0) +G2(J)*G3(J) )
331	AP=DENOM(G3(J)(G1(J)-1./UBAR0) +G2(J)*G4 )
332	C CPM1 AND CMM1 ARE THE CPLUS AND CMINUS TERMS EVALUATED
333	C AT THE TOP OF THE LAYER, THAT IS LOWER OPTICAL DEPTH TAU(J)
334	CPM1(J)=AP*TAU(J)
335	CMM1(J)=AM*TAU(J)
336	C CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
337	C BOTTOM OF THE LAYER. THAT IS AT HIGHER OPTICAL DEPTH TAU(J+1)
338	CP(J)=AP*TAU(J+1)
339	CM(J)=AM*TAU(J+1)
340	END DO
341
342	C NOW CALCULATE THE EXPONENTIAL TERMS NEEDED
343	C FOR THE TRIDIAGONAL ROTATED LAYERED METHOD
344	C WARNING IF DTAU(J) IS GREATER THAN ABOUT 35
345	C WE CLIPP IT TO AVOID OVERFLOW.
346	C EXP (TAU) - EXP(-TAU) WILL BE NONSENSE THIS IS
347	C CORRECTED IN THE DSOLVER ROUTINE. ??FLAG?
348	DO J=1,NAYER
349	c EXPTRM(J) = MIN(35.,LAMDA(J)*DTAU(J))
350	EP=EXP(MIN(35.0_8,LAMDA(J)*DTAU(J)))
351	EM=1.0/EP
352	AM = GAMA(J)*EM
353	E1(J)=EP+AM
354	E2(J)=EP-AM
355	AP = GAMA(J)*EP
356	E3(J)=AP+EM
357	E4(J)=AP-EM
358	END DO
359
360	CALL DSOLVER(NAYER,GAMA,CP,CM,CPM1,CMM1
361	& ,E1,E2,E3,E4,0.0_8,BSURF,RSF,XK1,XK2)
362
363	C EVALUATE THE NAYER FLUXES THROUGH THE NAYER LAYERS
364	C USE THE TOP (TAU=0) OPTICAL DEPTH EXPRESSIONS TO EVALUATE FP AND FM
365	C AT THE THE TOP OF EACH LAYER,J = LEVEL J
366	DO J=1,NAYER
367	FP(J)= XK1(J) + GAMA(J)*XK2(J) + CPM1(J)
368	FM(J)= GAMA(J)*XK1(J) + XK2(J) + CMM1(J)
369	END DO
370
371	C USE EXPRESSION FOR BOTTOM FLUX TO GET THE FP AND FM AT NLEV
372	c Optimization: no need this step since result of last
373	c loop at about EP above give this
374	c EP=EXP(EXPTRM(NAYER))
375	c EM=1.0/EP
376	FP(NLEV)=XK1(NAYER)EP + XK2(NAYER)AM + CP(NAYER)
377	FM(NLEV)=XK1(NAYER)AP + XK2(NAYER)EM + CM(NAYER)
378
379	C ADD THE DIRECT FLUX TERM TO THE DOWNWELLING RADIATION, LIOU 182
380	DO J=1,NLEV
381	FM(J)=FM(J)+UBAR0*TAU(J)
382	END DO
383
384	RETURN
385	END
386
387	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
388
389	SUBROUTINE DSOLVER(NL,GAMA,CP,CM,CPM1,CMM1,E1,E2,E3,E4,BTOP,
390	* BSURF,RSF,XK1,XK2)
391
392	C DOUBLE PRECISION VERSION OF SOLVER
393
394	cc PARAMETER (NMAX=201)
395	cc AS+JBM 03/2010
396	IMPLICIT REAL*8 (A-H,O-Z)
397	DIMENSION GAMA(NL),CP(NL),CM(NL),CPM1(NL),CMM1(NL),XK1(NL),
398	* XK2(NL),E1(NL),E2(NL),E3(NL),E4(NL)
399	cc AS+JBM 03/2010
400	cc DIMENSION AF(NMAX),BF(NMAX),CF(NMAX),DF(NMAX),XK(NMAX)
401	DIMENSION AF(2NL),BF(2NL),CF(2NL),DF(2NL),XK(2*NL)
402
403	C*********************************************************
404	C* THIS SUBROUTINE SOLVES FOR THE COEFFICIENTS OF THE *
405	C* TWO STREAM SOLUTION FOR GENERAL BOUNDARY CONDITIONS *
406	C* NO ASSUMPTION OF THE DEPENDENCE ON OPTICAL DEPTH OF *
407	C* C-PLUS OR C-MINUS HAS BEEN MADE. *
408	C* NL = NUMBER OF LAYERS IN THE MODEL *
409	C* CP = C-PLUS EVALUATED AT TAO=0 (TOP) *
410	C* CM = C-MINUS EVALUATED AT TAO=0 (TOP) *
411	C* CPM1 = C-PLUS EVALUATED AT TAOSTAR (BOTTOM) *
412	C* CMM1 = C-MINUS EVALUATED AT TAOSTAR (BOTTOM) *
413	C* EP = EXP(LAMDADTAU)
414	C* EM = 1/EP *
415	C* E1 = EP + GAMA EM
416	C* E2 = EP - GAMA EM
417	C* E3 = GAMAEP + EM
418	C* E4 = GAMAEP - EM
419	C* BTOP = THE DIFFUSE RADIATION INTO THE MODEL AT TOP *
420	C* BSURF = THE DIFFUSE RADIATION INTO THE MODEL AT *
421	C* THE BOTTOM: INCLUDES EMMISION AND REFLECTION *
422	C* OF THE UNATTENUATED PORTION OF THE DIRECT *
423	C* BEAM. BSTAR+RSFFOEXP(-TAOSTAR/U0) *
424	C* RSF = REFLECTIVITY OF THE SURFACE *
425	C* XK1 = COEFFICIENT OF THE POSITIVE EXP TERM *
426	C* XK2 = COEFFICIENT OF THE NEGATIVE EXP TERM *
427	C*********************************************************
428	C THIS ROUTINE CALLS ROUTINE DTRIDGL TO SOLVE TRIDIAGONAL
429	C SYSTEMS
430	C======================================================================C
431
432	L=2*NL
433
434	C **********MIXED COEFFICENTS********
435	C THIS VERSION AVOIDS SINGULARITIES ASSOC.
436	C WITH W0=0 BY SOLVING FOR XK1+XK2, AND XK1-XK2.
437
438	AF(1) = 0.0
439	BF(1) = GAMA(1)+1.
440	CF(1) = GAMA(1)-1.
441	DF(1) = BTOP-CMM1(1)
442	N = 0
443	LM2 = L-2
444
445	C EVEN TERMS
446
447	DO I=2,LM2,2
448	N = N+1
449	AF(I) = (E1(N)+E3(N))*(GAMA(N+1)-1.)
450	BF(I) = (E2(N)+E4(N))*(GAMA(N+1)-1.)
451	CF(I) = 2.0(1.-GAMA(N+1)*2)
452	DF(I) = (GAMA(N+1)-1.) * (CPM1(N+1) - CP(N)) +
453	* (1.-GAMA(N+1))* (CM(N)-CMM1(N+1))
454	END DO
455
456	N = 0
457	LM1 = L-1
458	DO I=3,LM1,2
459	N = N+1
460	AF(I) = 2.0(1.-GAMA(N)*2)
461	BF(I) = (E1(N)-E3(N))*(1.+GAMA(N+1))
462	CF(I) = (E1(N)+E3(N))*(GAMA(N+1)-1.)
463	DF(I) = E3(N)(CPM1(N+1) - CP(N)) + E1(N)(CM(N) - CMM1(N+1))
464	END DO
465
466	AF(L) = E1(NL)-RSF*E3(NL)
467	BF(L) = E2(NL)-RSF*E4(NL)
468	CF(L) = 0.0
469	DF(L) = BSURF-CP(NL)+RSF*CM(NL)
470
471	CALL DTRIDGL(L,AF,BF,CF,DF,XK)
472
473	C *UNMIX THE COEFFICIENTS**
474
475	DO 28 N=1,NL
476	XK1(N) = XK(2N-1)+XK(2N)
477	XK2(N) = XK(2N-1)-XK(2N)
478
479	C NOW TEST TO SEE IF XK2 IS REALLY ZERO TO THE LIMIT OF THE
480	C MACHINE ACCURACY = 1 .E -30
481	C XK2 IS THE COEFFICIENT OF THE GROWING EXPONENTIAL AND MUST
482	C BE TREATED CAREFULLY
483
484	IF(XK2(N) .EQ. 0.0) GO TO 28
485	IF (ABS (XK2(N)/XK(2*N-1)) .LT. 1.E-30) XK2(N)=0.0
486
487	28 CONTINUE
488
489	RETURN
490	END
491
492	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
493
494	SUBROUTINE DTRIDGL(L,AF,BF,CF,DF,XK)
495
496	C DOUBLE PRECISION VERSION OF TRIDGL
497
498	cc AS+JBM 03/2010 : OBSOLETE MAINTENANT
499	cc PARAMETER (NMAX=201)
500	IMPLICIT REAL*8 (A-H,O-Z)
501	DIMENSION AF(L),BF(L),CF(L),DF(L),XK(L)
502	cc AS+JBM 03/2010 : OBSOLETE MAINTENANT
503	cc DIMENSION AS(NMAX),DS(NMAX)
504	DIMENSION AS(L),DS(L)
505
506	C* THIS SUBROUTINE SOLVES A SYSTEM OF TRIDIAGIONAL MATRIX
507	C* EQUATIONS. THE FORM OF THE EQUATIONS ARE:
508	C* A(I)X(I-1) + B(I)X(I) + C(I)*X(I+1) = D(I)
509	C* WHERE I=1,L LESS THAN 103.
510	C* ..............REVIEWED -CP........
511
512	C======================================================================C
513
514	AS(L) = AF(L)/BF(L)
515	DS(L) = DF(L)/BF(L)
516
517	DO I=2,L
518	X = 1./(BF(L+1-I) - CF(L+1-I)*AS(L+2-I))
519	AS(L+1-I) = AF(L+1-I)*X
520	DS(L+1-I) = (DF(L+1-I)-CF(L+1-I)DS(L+2-I))X
521	END DO
522
523	XK(1)=DS(1)
524	DO I=2,L
525	XK(I) = DS(I)-AS(I)*XK(I-1)
526	END DO
527
528	RETURN
529	END
530

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format