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

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
[1047]	6	use dimradmars_mod, only: sunfr, ndlo2, nsun, ndlon, nflev
	7	use yomlw_h, only: nlaylte
[38]	8	IMPLICIT NONE
	9	C
[1047]	10	!#include "dimensions.h"
	11	!#include "dimphys.h"
	12	!#include "dimradmars.h"
[38]	13	#include "callkeys.h"
[1047]	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"
[38]	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
[1047]	42	c in dimradmars_mod , in each layer, for one of
[38]	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 )
[1047]	383	use dimradmars_mod, only: ndlo2
[38]	384	implicit none
	385	C
[1047]	386	!#include "dimensions.h"
	387	!#include "dimphys.h"
	388	!#include "dimradmars.h"
[38]	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