Context Navigation

swr_fouquart.F @ 1242

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

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

Download in other formats:

Original Format