Context Navigation

swr_fouquart.F @ 2156

Last change on this file since 2156 was 1266, checked in by aslmd, 11 years ago

LMDZ.MARS
IMPORTANT CHANGE

Remove all reference/use of nlayermx and dimphys.h
Made use of automatic arrays whenever arrays are needed with dimension nlayer
Remove lots of obsolete reference to dimensions.h
Converted iono.h and param_v4.h into corresponding modules

(with embedded subroutine to allocate arrays)
(no arrays allocated if thermosphere not used)

Deleted param.h and put contents into module param_v4_h
Adapted testphys1d, newstart, etc...
Made DATA arrays in param_read to be initialized by subroutine

fill_data_thermos in module param_v4_h

Optimized computations in paramfoto_compact (twice less dlog10 calculations)
Checked consistency before/after modification in debug mode
Checked performance is not impacted (same as before)

File size: 16.4 KB

Line
1	SUBROUTINE SWR_FOUQUART ( 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,
7	& ndlon, nflev, naerkind
8	use yomlw_h, only: nlaylte
9	IMPLICIT NONE
10	C
11	#include "callkeys.h"
12
13	C
14	C SWR - Continuum scattering computations
15	C
16	C PURPOSE.
17	C --------
18	C Computes the reflectivity and transmissivity in case oF
19	C Continuum scattering
20	c F. Forget (1999)
21	c
22	c BASED ON MORCRETTE EARTH MODEL
23	C (See radiation's part of the ecmwf research department
24	C documentation, and Fouquart and BonneL (1980)
25	C
26	C IMPLICIT ARGUMENTS :
27	C --------------------
28	C
29	C ==== INPUTS ===
30	c
31	c KDLON : number of horizontal grid points
32	c KFLEV : number of vertical layers
33	c KNU : Solar band # (1 or 2)
34	c aerosol aerosol extinction optical depth
35	c at reference wavelength "longrefvis" set
36	c in dimradmars_mod , in each layer, for one of
37	c the "naerkind" kind of aerosol optical properties.
38	c albedo hemispheric surface albedo
39	c albedo (i,1) : mean albedo for solar band#1
40	c (see below)
41	c albedo (i,2) : mean albedo for solar band#2
42	c (see below)
43	c PDSIG layer thickness in sigma coordinates
44	c PPSOL Surface pressure (Pa)
45	c PRMU: cos of solar zenith angle (=1 when sun at zenith)
46	c (CORRECTED for high zenith angle (atmosphere), unlike mu0)
47	c PSEC =1./PRMU
48
49	C ==== OUTPUTS ===
50	c
51	c PFD : downward flux in spectral band #INU in a given mesh
52	c (normalized to the total incident flux at the top of the atmosphere)
53	c PFU : upward flux in specatral band #INU in a given mesh
54	c (normalized to the total incident flux at the top of the atmosphere)
55	C
56	C
57	C METHOD.
58	C -------
59	C
60	C Computes continuum fluxes corresponding to aerosoL
61	C Or/and rayleigh scattering (no molecular gas absorption)
62	C
63	C-----------------------------------------------------------------------
64	C
65	C
66	C-----------------------------------------------------------------------
67	C
68
69	C ARGUMENTS
70	C ---------
71	INTEGER KDLON, KFLEV, KNU
72	REAL aerosol(NDLO2,KFLEV,naerkind), albedo(NDLO2,2),
73	S PDSIG(NDLO2,KFLEV),PSEC(NDLO2)
74
75	REAL QVISsQREF3d(NDLO2,KFLEV,nsun,naerkind)
76	REAL omegaVIS3d(NDLO2,KFLEV,nsun,naerkind)
77	REAL gVIS3d(NDLO2,KFLEV,nsun,naerkind)
78
79	REAL PPSOL(NDLO2)
80	REAL PFD(NDLO2,KFLEV+1),PFU(NDLO2,KFLEV+1)
81	REAL PRMU(NDLO2)
82
83	C LOCAL ARRAYS
84	C ------------
85
86	INTEGER jk,ja,jl,jae, jkl,jklp1,jkm1,jaj
87	REAL ZTRAY, ZRATIO,ZGAR, ZFF
88	real zfacoa,zcorae
89	real ZMUE, zgap,zbmu0, zww,zto,zden,zmu1,zbmu1,zden1,zre11
90
91	REAL ZC1I(NDLON,NFLEV+1), ZGG(NDLON), ZREF(NDLON)
92	S , ZRE1(NDLON), ZRE2(NDLON)
93	S , ZRMUZ(NDLON), ZRNEB(NDLON), ZR21(NDLON)
94	S , ZR23(NDLON), ZSS1(NDLON), ZTO1(NDLON), ZTR(NDLON,2,NFLEV+1)
95	S , ZTR1(NDLON), ZTR2(NDLON), ZW(NDLON)
96
97	REAL ZRAY1(NDLO2,NFLEV+1), ZRAY2(NDLO2,NFLEV+1)
98	s , ZREFZ(NDLO2,2,NFLEV+1)
99	S , ZRMUE(NDLO2,NFLEV+1)
100	S , ZCGAZ(NDLO2,NFLEV),ZPIZAZ(NDLO2,NFLEV),ZTAUAZ(NDLO2,NFLEV)
101
102	REAL ZRAYL(NDLON)
103	S , ZRJ(NDLON,6,NFLEV+1)
104	S , ZRK(NDLON,6,NFLEV+1)
105	S , ZTRA1(NDLON,NFLEV+1), ZTRA2(NDLON,NFLEV+1)
106
107	c Function
108	c --------
109	real CVMGT
110
111	C --------------------------------
112	C OPTICAL PARAMETERS FOR AEROSOLS
113	C -------------------------------
114	C
115	DO JK = 1 , nlaylte+1
116	DO JA = 1 , 6
117	DO JL = 1 , KDLON
118	ZRJ(JL,JA,JK) = 0.
119	ZRK(JL,JA,JK) = 0.
120	END DO
121	END DO
122	END DO
123
124	c Computing TOTAL single scattering parameters by adding
125	c properties of all the NAERKIND kind of aerosols
126
127	DO JK = 1 , nlaylte
128	DO JL = 1 , KDLON
129	ZCGAZ(JL,JK) = 0.
130	ZPIZAZ(JL,JK) = 0.
131	ZTAUAZ(JL,JK) = 0.
132	END DO
133	DO 106 JAE=1,naerkind
134	DO 105 JL = 1 , KDLON
135	c Mean Extinction optical depth in the spectral band
136	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
137	ZTAUAZ(JL,JK)=ZTAUAZ(JL,JK)
138	S +aerosol(JL,JK,JAE)*QVISsQREF3d(JL,JK,KNU,JAE)
139	c Single scattering albedo
140	c ~~~~~~~~~~~~~~~~~~~~~~~~
141	ZPIZAZ(JL,JK)=ZPIZAZ(JL,JK)+aerosol(JL,JK,JAE)*
142	S QVISsQREF3d(JL,JK,KNU,JAE)*
143	& omegaVIS3d(JL,JK,KNU,JAE)
144	c Assymetry factor
145	c ~~~~~~~~~~~~~~~~
146	ZCGAZ(JL,JK) = ZCGAZ(JL,JK) +aerosol(JL,JK,JAE)*
147	S QVISsQREF3d(JL,JK,KNU,JAE)*
148	& omegaVIS3d(JL,JK,KNU,JAE)*gVIS3d(JL,JK,KNU,JAE)
149	105 CONTINUE
150	106 CONTINUE
151	END DO
152	C
153	DO JK = 1 , nlaylte
154	DO JL = 1 , KDLON
155	ZCGAZ(JL,JK) = CVMGT( 0., ZCGAZ(JL,JK) / ZPIZAZ(JL,JK),
156	S (ZPIZAZ(JL,JK).EQ.0) )
157	ZPIZAZ(JL,JK) = CVMGT( 1., ZPIZAZ(JL,JK) / ZTAUAZ(JL,JK),
158	S (ZTAUAZ(JL,JK).EQ.0) )
159	END DO
160	END DO
161
162	C --------------------------------
163	C INCLUDING RAYLEIGH SCATERRING
164	C -------------------------------
165	if (rayleigh) then
166
167	call swrayleigh(kdlon,knu,ppsol,prmu,ZRAYL)
168
169	c Modifying mean aerosol parameters to account rayleigh scat by gas:
170
171	DO JK = 1 , nlaylte
172	DO JL = 1 , KDLON
173	c Rayleigh opacity in each layer :
174	ZTRAY = ZRAYL(JL) * PDSIG(JL,JK)
175	c ratio Tau(rayleigh) / Tau (total)
176	ZRATIO = ZTRAY / (ZTRAY + ZTAUAZ(JL,JK))
177	ZGAR = ZCGAZ(JL,JK)
178	ZFF = ZGAR * ZGAR
179	ZTAUAZ(JL,JK)=ZTRAY+ZTAUAZ(JL,JK)(1.-ZPIZAZ(JL,JK)ZFF)
180	ZCGAZ(JL,JK) = ZGAR * (1. - ZRATIO) / (1. + ZGAR)
181	ZPIZAZ(JL,JK) =ZRATIO+(1.-ZRATIO)ZPIZAZ(JL,JK)(1.-ZFF)
182	S / (1. -ZPIZAZ(JL,JK) * ZFF)
183	END DO
184	END DO
185	end if
186
187
188	C ----------------------------------------------
189	C TOTAL EFFECTIVE CLOUDINESS ABOVE A GIVEN LEVEL
190	C ----------------------------------------------
191	C
192	200 CONTINUE
193
194	DO JL = 1 , KDLON
195	ZR23(JL) = 0.
196	ZC1I(JL,nlaylte+1) = 0.
197	END DO
198
199	DO JK = 1 , nlaylte
200	JKL = nlaylte+1 - JK
201	JKLP1 = JKL + 1
202	DO JL = 1 , KDLON
203	ZFACOA = 1.-ZPIZAZ(JL,JKL)ZCGAZ(JL,JKL)ZCGAZ(JL,JKL)
204	ZCORAE = ZFACOA * ZTAUAZ(JL,JKL) * PSEC(JL)
205	ZR21(JL) = EXP(-ZCORAE )
206	ZSS1(JL) = 1.0-ZR21(JL)
207	ZC1I(JL,JKL) = 1.0-(1.0-ZSS1(JL))*(1.0-ZC1I(JL,JKLP1))
208	END DO
209	END DO
210
211	C -----------------------------------------------
212	C REFLECTIVITY/TRANSMISSIVITY FOR PURE SCATTERING
213	C -----------------------------------------------
214	C
215	DO JL = 1 , KDLON
216	ZRAY1(JL,nlaylte+1) = 0.
217	ZRAY2(JL,nlaylte+1) = 0.
218	ZREFZ(JL,2,1) = albedo(JL,KNU)
219	ZREFZ(JL,1,1) = albedo(JL,KNU)
220	ZTRA1(JL,nlaylte+1) = 1.
221	ZTRA2(JL,nlaylte+1) = 1.
222	END DO
223
224	DO JK = 2 , nlaylte+1
225	JKM1 = JK-1
226	DO 342 JL = 1 , KDLON
227	ZRNEB(JL)= 1.e-5 ! used to be "cloudiness" (PCLDSW in Morcrette)
228
229	ZRE1(JL)=0.
230	ZTR1(JL)=0.
231	ZRE2(JL)=0.
232	ZTR2(JL)=0.
233
234	C EQUIVALENT ZENITH ANGLE
235	c ~~~~~~~~~~~~~~~~~~~~~~~
236	ZMUE = (1.-ZC1I(JL,JK)) * PSEC(JL)
237	S + ZC1I(JL,JK) * 1.66
238	ZRMUE(JL,JK) = 1./ZMUE
239
240	C ------------------------------------------------------------------
241	C REFLECT./TRANSMISSIVITY DUE TO AEROSOLS (and rayleigh ?)
242	C ------------------------------------------------------------------
243
244	ZGAP = ZCGAZ(JL,JKM1)
245	ZBMU0 = 0.5 - 0.75 * ZGAP / ZMUE
246	ZWW =ZPIZAZ(JL,JKM1)
247	ZTO = ZTAUAZ(JL,JKM1)
248	ZDEN = 1. + (1. - ZWW + ZBMU0 * ZWW) * ZTO * ZMUE
249	S + (1-ZWW) * (1. - ZWW +2.ZBMU0ZWW)ZTOZTOZMUEZMUE
250	ZRAY1(JL,JKM1) = ZBMU0 * ZWW * ZTO * ZMUE / ZDEN
251	ZTRA1(JL,JKM1) = 1. / ZDEN
252	C
253	ZMU1 = 0.5
254	ZBMU1 = 0.5 - 0.75 * ZGAP * ZMU1
255	ZDEN1= 1. + (1. - ZWW + ZBMU1 * ZWW) * ZTO / ZMU1
256	S + (1-ZWW) * (1. - ZWW +2.ZBMU1ZWW)ZTOZTO/ZMU1/ZMU1
257	ZRAY2(JL,JKM1) = ZBMU1 * ZWW * ZTO / ZMU1 / ZDEN1
258	ZTRA2(JL,JKM1) = 1. / ZDEN1
259
260	ZGG(JL) = ZCGAZ(JL,JKM1)
261	ZW(JL) =ZPIZAZ(JL,JKM1)
262	ZREF(JL) = ZREFZ(JL,1,JKM1)
263	ZRMUZ(JL) = ZRMUE(JL,JK)
264	ZTO1(JL) = ZTAUAZ(JL,JKM1)/ZPIZAZ(JL,JKM1)
265
266	342 CONTINUE
267
268	C
269	CALL DEDD ( KDLON
270	S , ZGG,ZREF,ZRMUZ,ZTO1,ZW
271	S , ZRE1,ZRE2,ZTR1,ZTR2 )
272	C
273	DO JL = 1 , KDLON
274	C
275	ZREFZ(JL,1,JK) = (1.-ZRNEB(JL)) * (ZRAY1(JL,JKM1)
276	S + ZREFZ(JL,1,JKM1) * ZTRA1(JL,JKM1)
277	S * ZTRA2(JL,JKM1)
278	S / (1.-ZRAY2(JL,JKM1)*ZREFZ(JL,1,JKM1)))
279	S + ZRNEB(JL) * ZRE2(JL)
280	C
281	ZTR(JL,1,JKM1) = ZRNEB(JL) * ZTR2(JL) + (ZTRA1(JL,JKM1)
282	S / (1.-ZRAY2(JL,JKM1)*ZREFZ(JL,1,JKM1)))
283	S * (1.-ZRNEB(JL))
284	C
285	ZREFZ(JL,2,JK) = (1.-ZRNEB(JL)) * (ZRAY1(JL,JKM1)
286	S + ZREFZ(JL,2,JKM1) * ZTRA1(JL,JKM1)
287	S * ZTRA2(JL,JKM1) )
288	S + ZRNEB(JL) * ZRE1(JL)
289	C
290	ZTR(JL,2,JKM1) = ZRNEB(JL) * ZTR1(JL)
291	S + ZTRA1(JL,JKM1) * (1.-ZRNEB(JL))
292	C
293	END DO
294	END DO
295	C
296	C
297	C ------------------------------------------------------------------
298	C
299	C * 3.5 REFLECT./TRANSMISSIVITY BETWEEN SURFACE AND LEVEL
300	C -------------------------------------------------
301	C
302	350 CONTINUE
303	C
304	IF (KNU.EQ.1) THEN
305	JAJ = 2
306	DO 351 JL = 1 , KDLON
307	ZRJ(JL,JAJ,nlaylte+1) = 1.
308	ZRK(JL,JAJ,nlaylte+1) = ZREFZ(JL, 1,nlaylte+1)
309	351 CONTINUE
310	C
311	DO 353 JK = 1 , nlaylte
312	JKL = nlaylte+1 - JK
313	JKLP1 = JKL + 1
314	DO 352 JL = 1 , KDLON
315	ZRE11= ZRJ(JL,JAJ,JKLP1) * ZTR(JL, 1,JKL)
316	ZRJ(JL,JAJ,JKL) = ZRE11
317	ZRK(JL,JAJ,JKL) = ZRE11 * ZREFZ(JL, 1,JKL)
318	352 CONTINUE
319	353 CONTINUE
320	354 CONTINUE
321	C
322	ELSE
323	C
324	DO 358 JAJ = 1 , 2
325	DO 355 JL = 1 , KDLON
326	ZRJ(JL,JAJ,nlaylte+1) = 1.
327	ZRK(JL,JAJ,nlaylte+1) = ZREFZ(JL,JAJ,nlaylte+1)
328	355 CONTINUE
329	C
330	DO 357 JK = 1 , nlaylte
331	JKL = nlaylte+1 - JK
332	JKLP1 = JKL + 1
333	DO 356 JL = 1 , KDLON
334	ZRE11= ZRJ(JL,JAJ,JKLP1) * ZTR(JL,JAJ,JKL)
335	ZRJ(JL,JAJ,JKL) = ZRE11
336	ZRK(JL,JAJ,JKL) = ZRE11 * ZREFZ(JL,JAJ,JKL)
337	356 CONTINUE
338	357 CONTINUE
339	358 CONTINUE
340	END IF
341
342	C
343	C
344	C
345	C ------------------------------------------------------------------
346	C ---------------
347	C DOWNWARD FLUXES
348	C ---------------
349	C
350	JAJ = 2
351
352	do JK = 1 , nlaylte+1
353	JKL = nlaylte+1 - JK + 1
354	DO JL = 1 , KDLON
355	PFD(JL,JKL) = ZRJ(JL,JAJ,JKL) * sunfr(KNU)
356	end do
357	end do
358	C
359	C -------------
360	C UPWARD FLUXES
361	C -------------
362	DO JK = 1 , nlaylte+1
363	DO JL = 1 , KDLON
364	c ZRK = upward flux / incident top flux
365	PFU(JL,JK) = ZRK(JL,JAJ,JK) * sunfr(KNU)
366	END DO
367	END DO
368
369	C
370	RETURN
371	END
372
373	CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
374
375	SUBROUTINE DEDD (KDLON,PGG,PREF,PRMUZ,PTO1,PW
376	S , PRE1,PRE2,PTR1,PTR2 )
377	use dimradmars_mod, only: ndlo2
378	implicit none
379	C
380	C
381	C**** DEDD - DELTA-EDDINGTON IN A CLOUDY LAYER
382	C
383	C PURPOSE.
384	C --------
385	C COMPUTES THE REFLECTIVITY AND TRANSMISSIVITY OF A CLOUDY
386	C LAYER USING THE DELTA-EDDINGTON'S APPROXIMATION.
387	C
388	C** INTERFACE.
389	C ----------
390	C DEDD IS CALLED BY SW.
391	C
392	C SUBROUTINE DEDD (KDLON,PGG,PREF,PRMUZ,PTO1,PW
393	C S , PRE1,PRE2,PTR1,PTR2 )
394	C
395	C EXPLICIT ARGUMENTS :
396	C --------------------
397	C PGG : (NDLON) ; ASSYMETRY FACTOR
398	C PREF : (NDLON) ; REFLECTIVITY OF THE UNDERLYING LAYER
399	C PRMUZ : (NDLON) ; COSINE OF SOLAR ZENITH ANGLE
400	C PTO1 : (NDLON) ; OPTICAL THICKNESS
401	C PW : (NDLON) ; SINGLE SCATTERING ALBEDO
402	C ==== OUTPUTS ===
403	C PRE1 : (NDLON) ; LAYER REFLECTIVITY ASSUMING NO
404	C ; REFLECTION FROM UNDERLYING LAYER
405	C PTR1 : (NDLON) ; LAYER TRANSMISSIVITY ASSUMING NO
406	C ; REFLECTION FROM UNDERLYING LAYER
407	C PRE2 : (NDLON) ; LAYER REFLECTIVITY ASSUMING
408	C ; REFLECTION FROM UNDERLYING LAYER
409	C PTR2 : (NDLON) ; LAYER TRANSMISSIVITY ASSUMING
410	C ; REFLECTION FROM UNDERLYING LAYER
411	C
412	C IMPLICIT ARGUMENTS : NONE
413	C --------------------
414	C
415	C METHOD.
416	C -------
417	C
418	C STANDARD DELTA-EDDINGTON LAYER CALCULATIONS.
419	C
420	C EXTERNALS.
421	C ----------
422	C
423	C NONE
424	C
425	C REFERENCE.
426	C ----------
427	C
428	C SEE RADIATION'S PART OF THE MODEL'S DOCUMENTATION AND
429	C ECMWF RESEARCH DEPARTMENT DOCUMENTATION OF THE "IN CORE MODEL"
430	C
431	C AUTHOR.
432	C -------
433	C JEAN-JACQUES MORCRETTE ECMWF
434	C
435	C MODIFICATIONS.
436	C --------------
437	C ORIGINAL : 88-12-15
438	C ------------------------------------------------------------------
439	C
440	C* 0.1 ARGUMENTS
441	C ---------
442	INTEGER KDLON
443	C
444	REAL PGG(NDLO2),PREF(NDLO2),PRMUZ(NDLO2),PTO1(NDLO2),PW(NDLO2)
445	REAL PRE1(NDLO2),PRE2(NDLO2),PTR1(NDLO2),PTR2(NDLO2)
446
447	c local
448	integer jl
449	real*8 ZFF,ZGP,ZTOP,ZWCP,ZDT,ZX1,ZWM,ZRM2,ZRK,ZX2,ZRP,ZALPHA
450	real*8 ZBETA,ZEXMU0,ZEXKP,ZEXKM,ZXP2P,ZXM2P,ZAP2B,ZAM2B
451	real*8 ZA11,ZA12,ZA13,ZA22,ZA21,ZA23,ZDENA,ZC1A,ZC2A
452	real*8 ZRI0A,ZRI1A,ZRI0B,ZRI1B
453	real*8 ZB21,ZB22,ZB23,ZDENB,ZC1B,ZC2B
454	real*8 ZRI0C,ZRI1C,ZRI0D,ZRI1D
455	C
456	C ------------------------------------------------------------------
457	C
458	C* 1. DELTA-EDDINGTON CALCULATIONS
459	C
460	100 CONTINUE
461	C
462	DO 131 JL = 1 , KDLON
463	C
464	C* 1.1 SET UP THE DELTA-MODIFIED PARAMETERS
465	C
466	110 CONTINUE
467	C
468	ZFF = PGG(JL)*PGG(JL)
469	ZGP = PGG(JL)/(1.+PGG(JL))
470	ZTOP = (1.- PW(JL) * ZFF) * PTO1(JL)
471	ZWCP = (1-ZFF)* PW(JL) /(1.- PW(JL) * ZFF)
472	ZDT = 2./3.
473	ZX1 = 1.-ZWCP*ZGP
474	ZWM = 1.-ZWCP
475	ZRM2 = PRMUZ(JL) * PRMUZ(JL)
476	ZRK = SQRT(3.ZWMZX1)
477	ZX2 = 4.(1.-ZRKZRK*ZRM2)
478	ZRP = SQRT(3.*ZWM/ZX1)
479	ZALPHA = 3.ZWCPZRM2(1.+ZGPZWM)/ZX2
480	ZBETA = 3.ZWCP PRMUZ(JL) (1.+3.ZGPZRM2ZWM)/ZX2
481	ZEXMU0 = EXP(-ZTOP/ PRMUZ(JL) )
482	ZEXKP = EXP(ZRK*ZTOP)
483	ZEXKM = 1./ZEXKP
484	ZXP2P = 1.+ZDT*ZRP
485	ZXM2P = 1.-ZDT*ZRP
486	ZAP2B = ZALPHA+ZDT*ZBETA
487	ZAM2B = ZALPHA-ZDT*ZBETA
488	C
489	C* 1.2 WITHOUT REFLECTION FROM THE UNDERLYING LAYER
490	C
491	120 CONTINUE
492	C
493	ZA11 = ZXP2P
494	ZA12 = ZXM2P
495	ZA13 = ZAP2B
496	ZA22 = ZXP2P*ZEXKP
497	ZA21 = ZXM2P*ZEXKM
498	ZA23 = ZAM2B*ZEXMU0
499	ZDENA = ZA11 * ZA22 - ZA21 * ZA12
500	ZC1A = (ZA22ZA13-ZA12ZA23)/ZDENA
501	ZC2A = (ZA11ZA23-ZA21ZA13)/ZDENA
502	ZRI0A = ZC1A+ZC2A-ZALPHA
503	ZRI1A = ZRP*(ZC1A-ZC2A)-ZBETA
504	PRE1(JL) = (ZRI0A-ZDT*ZRI1A)/ PRMUZ(JL)
505	ZRI0B = ZC1AZEXKM+ZC2AZEXKP-ZALPHA*ZEXMU0
506	ZRI1B = ZRP(ZC1AZEXKM-ZC2AZEXKP)-ZBETAZEXMU0
507	PTR1(JL) = ZEXMU0+(ZRI0B+ZDT*ZRI1B)/ PRMUZ(JL)
508	C
509	C* 1.3 WITH REFLECTION FROM THE UNDERLYING LAYER
510	C
511	130 CONTINUE
512	C
513	ZB21 = ZA21- PREF(JL) ZXP2PZEXKM
514	ZB22 = ZA22- PREF(JL) ZXM2PZEXKP
515	ZB23 = ZA23- PREF(JL) ZEXMU0(ZAP2B - PRMUZ(JL) )
516	ZDENB = ZA11 * ZB22 - ZB21 * ZA12
517	ZC1B = (ZB22ZA13-ZA12ZB23)/ZDENB
518	ZC2B = (ZA11ZB23-ZB21ZA13)/ZDENB
519	ZRI0C = ZC1B+ZC2B-ZALPHA
520	ZRI1C = ZRP*(ZC1B-ZC2B)-ZBETA
521	PRE2(JL) = (ZRI0C-ZDT*ZRI1C) / PRMUZ(JL)
522	ZRI0D = ZC1BZEXKM + ZC2BZEXKP - ZALPHA*ZEXMU0
523	ZRI1D = ZRP * (ZC1BZEXKM - ZC2BZEXKP) - ZBETA*ZEXMU0
524	PTR2(JL) = ZEXMU0 + (ZRI0D + ZDT*ZRI1D) / PRMUZ(JL)
525	C
526	131 CONTINUE
527	RETURN
528	END
529

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

Download in other formats:

Original Format