Context Navigation

swr_toon.F @ 1233

Last change on this file since 1233 was 1047, checked in by emillour, 11 years ago

Mars GCM:

IMPORTANT CHANGE: Removed all reference/use of ngridmx (dimphys.h) in routines (necessary prerequisite to using parallel dynamics); in most cases this just means adding 'ngrid' as routine argument, and making local saved variables allocatable (and allocated at first call). In the process, had to convert many *.h files to equivalent modules: yomaer.h => yomaer_h.F90 , surfdat.h => surfdat_h.F90 , comsaison.h => comsaison_h.F90 , yomlw.h => yomlw_h.F90 , comdiurn.h => comdiurn_h.F90 , dimradmars.h => dimradmars_mod.F90 , comgeomfi.h => comgeomfi_h.F90, comsoil.h => comsoil_h.F90 , slope.h => slope_mod.F90
Also updated EOF routines, everything is now in eofdump_mod.F90
Removed unused routine lectfux.F (in dyn3d)

File size: 17.7 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	use dimradmars_mod, only: sunfr, ndlo2, nsun, nflev, ndlon
7	use yomlw_h, only: nlaylte
8
9	IMPLICIT NONE
10	C
11	!#include "dimensions.h"
12	!#include "dimphys.h"
13	!#include "dimradmars.h"
14	#include "callkeys.h"
15	! naerkind is set in scatterers.h (built when compiling with makegcm -s #)
16	#include"scatterers.h"
17	!#include "yomaer.h"
18	!#include "yomlw.h"
19
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	RETURN
216	END
217
218	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
219
220	SUBROUTINE GFLUXV(NAYER,DTDEL,WDEL,CDEL,UBAR0,RSF
221	& ,FP,FM)
222	IMPLICIT NONE
223	C THIS SUBROUTINE TAKES THE OPTICAL CONSTANTS AND BOUNDARY CONDITONS
224	C FOR THE VISIBLE FLUX AT ONE WAVELENGTH AND SOLVES FOR THE FLUXES AT
225	C THE LEVELS. THIS VERSION IS SET UP TO WORK WITH LAYER OPTICAL DEPTHS
226	C MEASURED FROM THE TOP OF EACH LAYER. (DTAU) TOP OF EACH LAYER HAS
227	C OPTICAL DEPTH TAU(N). IN THIS SUB LEVEL N IS ABOVE LAYER N. THAT IS LAYER N
228	C HAS LEVEL N ON TOP AND LEVEL N+1 ON BOTTOM. OPTICAL DEPTH INCREASES
229	C FROM TOP TO BOTTOM. SEE C.P. MCKAY, TGM NOTES.
230	C THIS SUBROUTINE DIFFERS FROM ITS IR CONTERPART IN THAT HERE WE SOLVE FOR
231	C THE FLUXES DIRECTLY USING THE GENERALIZED NOTATION OF MEADOR AND WEAVOR
232	C J.A.S., 37, 630-642, 1980.
233	C THE TRI-DIAGONAL MATRIX SOLVER IS DSOLVER AND IS DOUBLE PRECISION SO MANY
234	C VARIABLES ARE PASSED AS SINGLE THEN BECOME DOUBLE IN DSOLVER
235
236	C THIS VERSION HAS BEEN MODIFIED BY TRAN THE TRUNG WITH:
237	C 1. Simplified input & output for swr.F subroutine in LMDZ.MARS gcm model
238	C 2. Use delta function to modify optical properties
239	C 3. Use delta-eddington G1, G2, G3 parameters
240	C 4. Optimized for speed
241
242	C INPUTS:
243	INTEGER NAYER
244	c NAYER = number of layer
245	c first layer is at top
246	c last layer is at bottom
247	REAL*8 DTDEL(NAYER), WDEL(NAYER), CDEL(NAYER)
248	c DTDEL = optical thickness of layer
249	c WDEL = single scattering of layer
250	c CDEL = assymetry parameter
251	REAL*8 UBAR0, RSF
252	c UBAR0 = absolute value of cosine of solar zenith angle
253	c RSF = surface albedo
254	C OUTPUTS:
255	REAL*8 FP(NAYER+1), FM(NAYER+1)
256	c FP = flux up
257	c FM = flux down
258	C PRIVATES:
259	INTEGER J,NL,NLEV
260	!!!! AS+JBM 03/2010 BUG BUG si trop niveaux verticaux (LES)
261	!!!! ET PAS BESOIN DE HARDWIRE SALE ICI !
262	!!!! CORRIGER CE BUG AMELIORE EFFICACITE ET FLEXIBILITE
263	!! PARAMETER (NL=201)
264	!! C THIS VALUE (201) MUST BE .GE. 2*NAYER
265	REAL*8 BSURF,AP,AM,DENOM,EM,EP,G4
266	!! REAL*8 W0(NL), COSBAR(NL), DTAU(NL), TAU(NL)
267	!! REAL*8 LAMDA(NL),XK1(NL),XK2(NL)
268	!! REAL*8 G1(NL),G2(NL),G3(NL)
269	!! REAL*8 GAMA(NL),CP(NL),CM(NL),CPM1(NL),CMM1(NL)
270	!! REAL*8 E1(NL),E2(NL),E3(NL),E4(NL)
271	REAL8 W0(2NAYER), COSBAR(2NAYER), DTAU(2NAYER), TAU(2*NAYER)
272	REAL8 LAMDA(2NAYER),XK1(2NAYER),XK2(2NAYER)
273	REAL8 G1(2NAYER),G2(2NAYER),G3(2NAYER)
274	REAL8 GAMA(2NAYER),CP(2NAYER),CM(2NAYER),CPM1(2*NAYER)
275	REAL8 CMM1(2NAYER)
276	REAL8 E1(2NAYER),E2(2NAYER),E3(2NAYER),E4(2*NAYER)
277
278	NL = 2*NAYER !!! AS+JBM 03/2010
279	NLEV = NAYER+1
280
281	C TURN ON THE DELTA-FUNCTION IF REQUIRED HERE
282	c TAU(1) = 0.0
283	c DO J=1,NAYER
284	c W0(J)=WDEL(J)
285	c COSBAR(J)=CDEL(J)
286	c DTAU(J)=DTDEL(J)
287	c TAU(J+1)=TAU(J)+DTAU(J)
288	c END DO
289	C FOR THE DELTA FUNCTION HERE...
290	TAU(1) = 0.0
291	DO J=1,NAYER
292	COSBAR(J)=CDEL(J)/(1.+CDEL(J))
293	W0(J)=1.-WDEL(J)CDEL(J)*2
294	DTAU(J)=DTDEL(J)*W0(J)
295	W0(J)=WDEL(J)(1.-CDEL(J)*2)/W0(J)
296	TAU(J+1)=TAU(J)+DTAU(J)
297	END DO
298
299	c Optimization, this is the major speed gain
300	TAU(1) = 1.0
301	do J = 2, NAYER+1
302	TAU(J) = EXP(-TAU(J)/UBAR0)
303	end do
304	BSURF = RSFUBAR0TAU(NLEV)
305	C WE GO WITH THE HEMISPHERIC CONSTANT APPROACH
306	C AS DEFINED BY M&W - THIS IS THE WAY THE IR IS DONE
307	DO J=1,NAYER
308	c Optimization: ALPHA not used
309	c ALPHA(J)=SQRT( (1.-W0(J))/(1.-W0(J)*COSBAR(J)) )
310	C SET OF CONSTANTS DETERMINED BY DOM
311	c G1(J)= (SQRT(3.)0.5)(2. - W0(J)*(1.+COSBAR(J)))
312	c G2(J)= (SQRT(3.)W0(J)0.5)*(1.-COSBAR(J))
313	c G3(J)=0.5(1.-SQRT(3.)COSBAR(J)*UBAR0)
314	c We use delta-Eddington instead
315	G1(J)=0.25(7.-W0(j)(4.+3*cosbar(j)))
316	G2(J)=-0.25(1.-W0(j)(4.-3*cosbar(j)))
317	G3(J)=0.5(1.-SQRT(3.)COSBAR(J)*UBAR0)
318	LAMDA(J)=SQRT(G1(J)2 - G2(J)2)
319	GAMA(J)=(G1(J)-LAMDA(J))/G2(J)
320	END DO
321
322	DO J=1,NAYER
323	G4=1.-G3(J)
324	DENOM=LAMDA(J)2 - 1./UBAR02
325	C NOTE THAT THE ALGORITHM DONOT ACCEPT UBAR0 .eq. 0
326	C THERE IS A POTENTIAL PROBLEM HERE IF W0=0 AND UBARV=UBAR0
327	C THEN DENOM WILL VANISH. THIS ONLY HAPPENS PHYSICALLY WHEN
328	C THE SCATTERING GOES TO ZERO
329	C PREVENT THIS WITH AN IF STATEMENT
330	IF ( DENOM .EQ. 0.) THEN
331	DENOM=1.E-6
332	END IF
333	DENOM = W0(J)/DENOM
334	AM=DENOM(G4 (G1(J)+1./UBAR0) +G2(J)*G3(J) )
335	AP=DENOM(G3(J)(G1(J)-1./UBAR0) +G2(J)*G4 )
336	C CPM1 AND CMM1 ARE THE CPLUS AND CMINUS TERMS EVALUATED
337	C AT THE TOP OF THE LAYER, THAT IS LOWER OPTICAL DEPTH TAU(J)
338	CPM1(J)=AP*TAU(J)
339	CMM1(J)=AM*TAU(J)
340	C CP AND CM ARE THE CPLUS AND CMINUS TERMS EVALUATED AT THE
341	C BOTTOM OF THE LAYER. THAT IS AT HIGHER OPTICAL DEPTH TAU(J+1)
342	CP(J)=AP*TAU(J+1)
343	CM(J)=AM*TAU(J+1)
344	END DO
345
346	C NOW CALCULATE THE EXPONENTIAL TERMS NEEDED
347	C FOR THE TRIDIAGONAL ROTATED LAYERED METHOD
348	C WARNING IF DTAU(J) IS GREATER THAN ABOUT 35
349	C WE CLIPP IT TO AVOID OVERFLOW.
350	C EXP (TAU) - EXP(-TAU) WILL BE NONSENSE THIS IS
351	C CORRECTED IN THE DSOLVER ROUTINE. ??FLAG?
352	DO J=1,NAYER
353	c EXPTRM(J) = MIN(35.,LAMDA(J)*DTAU(J))
354	EP=EXP(MIN(35.0_8,LAMDA(J)*DTAU(J)))
355	EM=1.0/EP
356	AM = GAMA(J)*EM
357	E1(J)=EP+AM
358	E2(J)=EP-AM
359	AP = GAMA(J)*EP
360	E3(J)=AP+EM
361	E4(J)=AP-EM
362	END DO
363
364	CALL DSOLVER(NAYER,GAMA,CP,CM,CPM1,CMM1
365	& ,E1,E2,E3,E4,0.0_8,BSURF,RSF,XK1,XK2)
366
367	C EVALUATE THE NAYER FLUXES THROUGH THE NAYER LAYERS
368	C USE THE TOP (TAU=0) OPTICAL DEPTH EXPRESSIONS TO EVALUATE FP AND FM
369	C AT THE THE TOP OF EACH LAYER,J = LEVEL J
370	DO J=1,NAYER
371	FP(J)= XK1(J) + GAMA(J)*XK2(J) + CPM1(J)
372	FM(J)= GAMA(J)*XK1(J) + XK2(J) + CMM1(J)
373	END DO
374
375	C USE EXPRESSION FOR BOTTOM FLUX TO GET THE FP AND FM AT NLEV
376	c Optimization: no need this step since result of last
377	c loop at about EP above give this
378	c EP=EXP(EXPTRM(NAYER))
379	c EM=1.0/EP
380	FP(NLEV)=XK1(NAYER)EP + XK2(NAYER)AM + CP(NAYER)
381	FM(NLEV)=XK1(NAYER)AP + XK2(NAYER)EM + CM(NAYER)
382
383	C ADD THE DIRECT FLUX TERM TO THE DOWNWELLING RADIATION, LIOU 182
384	DO J=1,NLEV
385	FM(J)=FM(J)+UBAR0*TAU(J)
386	END DO
387
388	RETURN
389	END
390
391	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
392
393	SUBROUTINE DSOLVER(NL,GAMA,CP,CM,CPM1,CMM1,E1,E2,E3,E4,BTOP,
394	* BSURF,RSF,XK1,XK2)
395
396	C DOUBLE PRECISION VERSION OF SOLVER
397
398	cc PARAMETER (NMAX=201)
399	cc AS+JBM 03/2010
400	IMPLICIT REAL*8 (A-H,O-Z)
401	DIMENSION GAMA(NL),CP(NL),CM(NL),CPM1(NL),CMM1(NL),XK1(NL),
402	* XK2(NL),E1(NL),E2(NL),E3(NL),E4(NL)
403	cc AS+JBM 03/2010
404	cc DIMENSION AF(NMAX),BF(NMAX),CF(NMAX),DF(NMAX),XK(NMAX)
405	DIMENSION AF(2NL),BF(2NL),CF(2NL),DF(2NL),XK(2*NL)
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	L=2*NL
437
438	C **********MIXED COEFFICENTS********
439	C THIS VERSION AVOIDS SINGULARITIES ASSOC.
440	C WITH W0=0 BY SOLVING FOR XK1+XK2, AND XK1-XK2.
441
442	AF(1) = 0.0
443	BF(1) = GAMA(1)+1.
444	CF(1) = GAMA(1)-1.
445	DF(1) = BTOP-CMM1(1)
446	N = 0
447	LM2 = L-2
448
449	C EVEN TERMS
450
451	DO I=2,LM2,2
452	N = N+1
453	AF(I) = (E1(N)+E3(N))*(GAMA(N+1)-1.)
454	BF(I) = (E2(N)+E4(N))*(GAMA(N+1)-1.)
455	CF(I) = 2.0(1.-GAMA(N+1)*2)
456	DF(I) = (GAMA(N+1)-1.) * (CPM1(N+1) - CP(N)) +
457	* (1.-GAMA(N+1))* (CM(N)-CMM1(N+1))
458	END DO
459
460	N = 0
461	LM1 = L-1
462	DO I=3,LM1,2
463	N = N+1
464	AF(I) = 2.0(1.-GAMA(N)*2)
465	BF(I) = (E1(N)-E3(N))*(1.+GAMA(N+1))
466	CF(I) = (E1(N)+E3(N))*(GAMA(N+1)-1.)
467	DF(I) = E3(N)(CPM1(N+1) - CP(N)) + E1(N)(CM(N) - CMM1(N+1))
468	END DO
469
470	AF(L) = E1(NL)-RSF*E3(NL)
471	BF(L) = E2(NL)-RSF*E4(NL)
472	CF(L) = 0.0
473	DF(L) = BSURF-CP(NL)+RSF*CM(NL)
474
475	CALL DTRIDGL(L,AF,BF,CF,DF,XK)
476
477	C *UNMIX THE COEFFICIENTS**
478
479	DO 28 N=1,NL
480	XK1(N) = XK(2N-1)+XK(2N)
481	XK2(N) = XK(2N-1)-XK(2N)
482
483	C NOW TEST TO SEE IF XK2 IS REALLY ZERO TO THE LIMIT OF THE
484	C MACHINE ACCURACY = 1 .E -30
485	C XK2 IS THE COEFFICIENT OF THE GROWING EXPONENTIAL AND MUST
486	C BE TREATED CAREFULLY
487
488	IF(XK2(N) .EQ. 0.0) GO TO 28
489	IF (ABS (XK2(N)/XK(2*N-1)) .LT. 1.E-30) XK2(N)=0.0
490
491	28 CONTINUE
492
493	RETURN
494	END
495
496	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
497
498	SUBROUTINE DTRIDGL(L,AF,BF,CF,DF,XK)
499
500	C DOUBLE PRECISION VERSION OF TRIDGL
501
502	cc AS+JBM 03/2010 : OBSOLETE MAINTENANT
503	cc PARAMETER (NMAX=201)
504	IMPLICIT REAL*8 (A-H,O-Z)
505	DIMENSION AF(L),BF(L),CF(L),DF(L),XK(L)
506	cc AS+JBM 03/2010 : OBSOLETE MAINTENANT
507	cc DIMENSION AS(NMAX),DS(NMAX)
508	DIMENSION AS(L),DS(L)
509
510	C* THIS SUBROUTINE SOLVES A SYSTEM OF TRIDIAGIONAL MATRIX
511	C* EQUATIONS. THE FORM OF THE EQUATIONS ARE:
512	C* A(I)X(I-1) + B(I)X(I) + C(I)*X(I+1) = D(I)
513	C* WHERE I=1,L LESS THAN 103.
514	C* ..............REVIEWED -CP........
515
516	C======================================================================C
517
518	AS(L) = AF(L)/BF(L)
519	DS(L) = DF(L)/BF(L)
520
521	DO I=2,L
522	X = 1./(BF(L+1-I) - CF(L+1-I)*AS(L+2-I))
523	AS(L+1-I) = AF(L+1-I)*X
524	DS(L+1-I) = (DF(L+1-I)-CF(L+1-I)DS(L+2-I))X
525	END DO
526
527	XK(1)=DS(1)
528	DO I=2,L
529	XK(I) = DS(I)-AS(I)*XK(I-1)
530	END DO
531
532	RETURN
533	END
534

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

Download in other formats:

Original Format