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

Rev	Line
[38]	1	SUBROUTINE SWR_FOUQUART ( KDLON, KFLEV, KNU
	2	S , aerosol,QVISsQREF3d,omegaVIS3d,gVIS3d
	3	& , albedo,PDSIG,PPSOL,PRMU,PSEC
	4	S , PFD,PFU )
	5
[1246]	6	use dimradmars_mod, only: sunfr, ndlo2, nsun,
	7	& ndlon, nflev, naerkind
[1047]	8	use yomlw_h, only: nlaylte
[38]	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
[1047]	36	c in dimradmars_mod , in each layer, for one of
[38]	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 )
[1047]	377	use dimradmars_mod, only: ndlo2
[38]	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