Context Navigation

swr_toon.F @ 3727

Last change on this file since 3727 was 3727, checked in by emillour, 2 months ago
Mars PCM: Code tidying: put routines in modules, remove useless "return" statements and remove obsolete and unused scopyi.F EM
File size: 17.9 KB

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

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

Download in other formats:

Original Format