Context Navigation

swr_toon.F @ 3807

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

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

Download in other formats:

Original Format