Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

sisvat_sno_albedo.F @ 4350

Last change on this file since 4350 was 3980, checked in by evignon, 3 years ago
je renomme la routine sisvat_sop.F pour faciliter les passerelles entre MAR-SISVAT et LMDZ-SISVAT
File size: 32.0 KB

Rev	Line
[3792]	1	subroutine SnOptP(jjtime)
	2
	3	C +------------------------------------------------------------------------+
	4	C \| MAR/SISVAT SnOptP 12-08-2019 MAR \|
	5	C \| SubRoutine SnOptP computes the Snow Pack optical Properties \|
	6	C +------------------------------------------------------------------------+
	7	C \| \|
	8	C \| PARAMETERS: knonv: Total Number of columns = \|
	9	C \| ^^^^^^^^^^ = Total Number of continental Grid Boxes \|
	10	C \| X Number of Mosaic Cell per Grid Box \|
	11	C \| \|
	12	C \| INPUT: isnoSV = total Nb of Ice/Snow Layers \|
	13	C \| ^^^^^ ispiSV = 0,...,nsno: Uppermost Superimposed Ice Layer \|
	14	C \| \|
	15	C \| \|
	16	C \| INPUT: G1snSV : Dendricity (<0) or Sphericity (>0) of Snow Layer \|
	17	C \| ^^^^^ G2snSV : Sphericity (>0) or Size of Snow Layer \|
	18	C \| agsnSV : Snow Age [day] \|
	19	C \| ro__SV : Snow/Soil Volumic Mass [kg/m3] \|
	20	C \| eta_SV : Water Content [m3/m3] \|
	21	C \| rusnSV : Surficial Water Thickness [kg/m2] .OR. [mm] \|
	22	C \| SWS_SV : Surficial Water Status \|
	23	C \| dzsnSV : Snow Layer Thickness [m] \|
	24	C \| \|
	25	C \| albssv : Soil Albedo [-] \|
	26	C \| zzsnsv : Snow Pack Thickness [m] \|
	27	C \| \|
	28	C \| OUTPUT: albisv : Snow/Ice/Water/Soil Integrated Albedo [-] \|
	29	C \| ^^^^^^ sEX_sv : Verticaly Integrated Extinction Coefficient \|
[3900]	30	C \| DOPsnSV : Snow Optical diameter [m] \|
[3792]	31	C \| \|
	32	C \| Internal Variables: \|
	33	C \| ^^^^^^^^^^^^^^^^^^ \|
	34	C \| SnOpSV : Snow Grain optical Size [m] \|
	35	C \| EX1_sv : Integrated Snow Extinction (0.3--0.8micr.m) \|
	36	C \| EX2_sv : Integrated Snow Extinction (0.8--1.5micr.m) \|
	37	C \| EX3_sv : Integrated Snow Extinction (1.5--2.8micr.m) \|
	38	C \| \|
	39	C \| METHODE: Calcul de la taille optique des grains ? partir de \|
	40	C \| ^^^^^^^ -leur type decrit par les deux variables descriptives \|
	41	C \| continues sur la plage -99/+99 passees en appel. \|
	42	C \| -la taille optique (1/10mm) des etoiles, \|
	43	C \| des grains fins et \|
	44	C \| des jeunes faces planes \|
	45	C \| \|
	46	C \| METHOD: Computation of the optical diameter of the grains \|
	47	C \| ^^^^^^ described with the CROCUS formalism G1snSV / G2snSV \|
	48	C \| \|
	49	C \| REFERENCE: Brun et al. 1989, J. Glaciol 35 pp. 333--342 \|
	50	C \| ^^^^^^^^^ Brun et al. 1992, J. Glaciol 38 pp. 13-- 22 \|
	51	C \| Eric Martin Sept.1996 \|
	52	C \| \|
	53	C \| \|
	54	C +------------------------------------------------------------------------+
	55
	56
	57
	58
	59	C +--Global Variables
	60	C + ================
	61
	62
	63	use VARphy
	64	use VAR_SV
	65	use VARdSV
	66	use VARxSV
	67	use VARySV
	68	use VARtSV
[3900]	69	USE surface_data, only: iflag_albcalc,correc_alb
[3792]	70
	71	IMPLICIT NONE
	72
	73
	74	C + -- INPUT
	75	integer jjtime
	76
	77	C +--Internal Variables
	78	C + ==================
	79
	80	real coalb1(knonv) ! weighted Coalbedo, Vis.
	81	real coalb2(knonv) ! weighted Coalbedo, nIR 1
	82	real coalb3(knonv) ! weighted Coalbedo, nIR 2
	83	real coalbm ! weighted Coalbedo, mean
	84	real sExt_1(knonv) ! Extinction Coeff., Vis.
	85	real sExt_2(knonv) ! Extinction Coeff., nIR 1
	86	real sExt_3(knonv) ! Extinction Coeff., nIR 2
	87	real SnOpSV(knonv, nsno) ! Snow Grain optical Size
	88	c #AG real agesno
	89
[3900]	90	integer isn ,ikl ,isn1, i
[3792]	91	real sbeta1,sbeta2,sbeta3,sbeta4,sbeta5
	92	real AgeMax,AlbMin,HSnoSV,HIceSV,doptmx,SignG1,Sph_OK
	93	real dalbed,dalbeS,dalbeW
[3900]	94	real bsegal,czemax,csegal,csza
[3792]	95	real RoFrez,DiffRo,SignRo,SnowOK,OpSqrt
	96	real albSn1,albIc1,a_SnI1,a_SII1
	97	real albSn2,albIc2,a_SnI2,a_SII2
	98	real albSn3,albIc3,a_SnI3,a_SII3
	99	real albSno,albIce,albSnI,albSII,albWIc,albmax
	100	real doptic,Snow_H,SIce_H,SnownH,SIcenH
	101	real exarg1,exarg2,exarg3,sign_0,sExt_0
	102	real albedo_old,albCor
	103	real ro_ave,dz_ave,minalb
[3900]	104	real l1min,l1max,l2min,l2max,l3min,l3max
	105	real l6min(6), l6max(6), albSn6(6), a_SII6(6)
	106	real lmintmp,lmaxtmp,albtmp
[3792]	107
	108	C +--Local DATA
	109	C + ============
	110
	111	C +--For the computation of the solar irradiance extinction in snow
	112	C + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	113	data sbeta1/0.0192/,sbeta2/0.4000/,sbeta3/0.1098/
	114	data sbeta4/1.0000/
	115	data sbeta5/2.00e1/
	116
	117	C +--Snow Age Maximum (Taiga, e.g. Col de Porte)
	118	C + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	119	data AgeMax /60.0/
	120	C +... AgeMax: Snow Age Maximum [day]
	121
	122	data AlbMin /0.94/
	123	C +... AlbMin: Albedo Minimum / visible (0.3--0.8 micrometers)
	124
	125	data HSnoSV /0.01/
	126	C +... HSnoSV: Snow Thickness through witch
	127	C + Albedo is interpolated to Ice Albedo
	128	data HIceSV /0.10/
	129	C +... HIceSV: Snow/Ice Thickness through witch
	130	C + Albedo is interpolated to Soil Albedo
	131	data doptmx /2.3e-3/
	132	C +... doptmx: Maximum optical Diameter (pi * R**2) [m]
	133	C +
	134	data czeMAX /0.173648178/ ! 80.deg (Segal et al., 1991 JAS)
	135
	136	data bsegal /4.00 / !
	137	data albmax /0.99 / ! Albedo max
	138
[3900]	139	C +-- wavelength limits [m] for each broad band
[3792]	140
[3900]	141	data l1min/400.0e-9/,l1max/800.0e-9/
	142	data l2min/800.0e-9/,l2max/1500.0e-9/
	143	data l3min/1500.0e-9/,l3max/2800.0e-9/
	144
	145	data l6min/185.0e-9,250.0e-9,400.0e-9,
	146	. 690.0e-9,1190.0e-9,2380.0e-9/
	147	data l6max/250.0e-9,400.0e-9,
	148	. 690.0e-9,1190.0e-9,2380.0e-9,4000.0e-9/
[3792]	149
	150
	151	C +--Snow Grain optical Size
	152	C + =======================
	153
	154	DO ikl=1,knonv
	155	DO isn=1,max(1,isnoSV(ikl))
	156
	157	G2snSV(ikl,isn) = max(epsi,G2snSV(ikl,isn))
	158	C +... Avoid non physical Values
	159
	160	SignG1 = sign(unun,G1snSV(ikl,isn))
	161	Sph_OK = max(zero,SignG1)
	162
	163	SnOpSV(ikl,isn) = 1.e-4 *
	164	C +... SI: (from 1/10 mm to m)
	165
	166
	167	C +--Contribution of Non Dendritic Snow
	168	C + ----------------------------------
	169
	170	. ( Sph_OK ( G2snSV(ikl,isn)G1snSV(ikl,isn)/G1_dSV
	171	. +max(demi*G2snSV(ikl,isn),DFcdSV)
	172	. *(unun-G1snSV(ikl,isn) /G1_dSV))
	173
	174
	175	C +--Contribution of Dendritic Snow
	176	C + ----------------------------------
	177
	178	. +(1.-Sph_OK)( -G1snSV(ikl,isn)DDcdSV /G1_dSV
	179	. +(unun+G1snSV(ikl,isn) /G1_dSV)
	180	. * (G2snSV(ikl,isn)*DScdSV /G1_dSV
	181	. +(unun-G2snSV(ikl,isn) /G1_dSV)
	182	. *DFcdSV )))
	183	SnOpSV(ikl,isn) = max(zero,SnOpSV(ikl,isn))
[3900]	184
	185	C + --For outputs (Etienne)
	186	C + ------------------------
	187	DOPsnSV(ikl,isn)=SnOpSV(ikl,isn)
[3792]	188	END DO
	189	END DO
	190
[3900]	191
	192
[3792]	193
	194	C +--Snow/Ice Albedo
	195	C + ===============
	196
	197
	198
	199	C +--Uppermost effective Snow Layer
	200	C + ------------------------------
	201
	202	DO ikl=1,knonv
	203
[3900]	204	isn = max(iun,isnoSV(ikl))
[3792]	205
	206	SignRo = sign(unun, rocdSV - ro__SV(ikl,isn))
	207	SnowOK = max(zero,SignRo) ! Ice Density Threshold
	208
	209	OpSqrt = sqrt(SnOpSV(ikl,isn))
	210
	211	cCA +--Correction of snow albedo for Antarctica/Greenland
	212	cCA --------------------------------------------------
[3900]	213
	214
	215	albCor = correc_alb
[3792]	216	c #GL albCor = 1.01
[3900]	217	c #AC albCor = 1.01
	218
	219
	220	IF (iflag_albcalc .GE. 1) THEN ! Albedo calculation according to Kokhanovsky and Zege 2004
	221
	222	dalbed = 0.0
	223	doptic=SnOpSV(ikl,isn)
	224	csza=coszSV(ikl)
	225
	226	CALL albedo_kokhanovsky(l1min,l1max,csza,doptic,albSn1)
	227	CALL albedo_kokhanovsky(l2min,l2max,csza,doptic,albSn2)
	228	CALL albedo_kokhanovsky(l3min,l3max,csza,doptic,albSn3)
	229
	230	DO i=1,6
	231	lmintmp=l6min(i)
	232	lmaxtmp=l6max(i)
	233	CALL albedo_kokhanovsky(lmintmp,lmaxtmp,csza,doptic,albtmp)
	234	albSn6(i)=albtmp
	235	ENDDO
	236
	237
	238	ELSE ! Default calculation in SISVAT
	239
	240	! Zenith Angle Correction (Segal et al., 1991, JAS 48, p.1025)
	241	!--------------------------- (Wiscombe & Warren, dec1980, JAS , p.2723)
	242	! (Warren, 1982, RG , p. 81)
	243	! -------------------------------------------------
	244
	245	dalbed = 0.0
	246
	247	csegal = max(czemax ,coszSV(ikl))
	248	c #cz dalbeS = ((bsegal+unun)/(unun+2.0bsegalcsegal)
	249	c #cz. - unun )*0.32
	250	c #cz. / bsegal
	251	c #cz dalbeS = max(dalbeS,zero)
	252	c #cz dalbed = dalbeS * min(1,isnoSV(ikl))
[3792]	253
[3900]	254	dalbeW =(0.64 - csegal )*0.0625 ! Warren 1982, RevGeo, fig.12b
	255	! 0.0625 = 5% * 1/0.8, p.81
	256	! 0.64 = cos(50)
	257	dalbed = dalbeW * min(1,isnoSV(ikl))
	258	!-------------------------------------------------------------------------
	259
[3792]	260	albSn1 = 0.96-1.580*OpSqrt
	261	albSn1 = max(albSn1,AlbMin)
	262
	263	albSn1 = max(albSn1,zero)
	264	albSn1 = min(albSn1*albCor,unun)
	265
	266	albSn2 = 0.95-15.40*OpSqrt
	267	albSn2 = max(albSn2,zero)
	268	albSn2 = min(albSn2*albCor,unun)
	269
	270	doptic = min(SnOpSV(ikl,isn),doptmx)
	271	albSn3 = 346.3doptic -32.31OpSqrt +0.88
	272	albSn3 = max(albSn3,zero)
	273	albSn3 = min(albSn3*albCor,unun)
	274
[3900]	275	albSn6(1:3)=albSn1
	276	albSn6(4:6)=albSn2
	277
[3792]	278	!snow albedo corection if wetsnow
	279	c #GL albSn1 = albSn1max(0.9,(1.-1.5eta_SV(ikl,isn)))
	280	c #GL albSn2 = albSn2max(0.9,(1.-1.5eta_SV(ikl,isn)))
	281	c #GL albSn3 = albSn3max(0.9,(1.-1.5eta_SV(ikl,isn)))
[3900]	282
	283	ENDIF
	284
[3792]	285
	286	albSno = So1dSV*albSn1
	287	. + So2dSV*albSn2
	288	. + So3dSV*albSn3
	289
	290	cXF
	291	minalb = (aI2dSV + (aI3dSV -aI2dSV)
	292	. * (ro__SV(ikl,isn)-ro_Ice)/(roSdSV-ro_Ice))
	293	minalb = min(aI3dSV,max(aI2dSV,minalb)) ! pure/firn albedo
	294
	295	SnowOK = SnowOK*max(zero,sign(unun, albSno-minalb))
	296	albSn1 = SnowOKalbSn1+(1.0-SnowOK)max(albSno,minalb)
	297	albSn2 = SnowOKalbSn2+(1.0-SnowOK)max(albSno,minalb)
	298	albSn3 = SnowOKalbSn3+(1.0-SnowOK)max(albSno,minalb)
[3900]	299	albSn6(:) = SnowOKalbSn6(:)+(1.0-SnowOK)max(albSno,minalb)
	300
[3792]	301
	302	c + ro < roSdSV \| min al > aI3dSV
	303	c + roSdSV < ro < rocdSV \| aI2dSV < min al < aI3dSV (fct of density)
	304
	305
	306	C +--Snow/Ice Pack Thickness
	307	C + -----------------------
	308
	309	isn = max(min(isnoSV(ikl) ,ispiSV(ikl)),0)
	310	Snow_H = zzsnsv(ikl,isnoSV(ikl))-zzsnsv(ikl,isn)
	311	SIce_H = zzsnsv(ikl,isnoSV(ikl))
	312	SnownH = Snow_H / HSnoSV
	313	SnownH = min(unun, SnownH)
	314	SIcenH = SIce_H / (HIceSV)
	315	SIcenH = min(unun, SIcenH)
	316
	317	C + The value of SnownH is set to 1 in case of ice lenses above
	318	C + 1m of dry snow (ro<600kg/m3) for using CROCUS albedo
	319
	320	c ro_ave = 0.
	321	c dz_ave = 0.
	322	c SnowOK = 1.
	323	c do isn = isnoSV(ikl),1,-1
	324	c ro_ave = ro_ave + ro__SV(ikl,isn) * dzsnSV(ikl,isn) * SnowOK
	325	c dz_ave = dz_ave + dzsnSV(ikl,isn) * SnowOK
	326	c SnowOK = max(zero,sign(unun,1.-dz_ave))
	327	c enddo
	328
	329	c ro_ave = ro_ave / max(dz_ave,epsi)
	330	c SnowOK = max(zero,sign(unun,600.-ro_ave))
	331	c SnownH = SnowOK + SnownH * (1. - SnowOK)
	332
	333	C +--Integrated Snow/Ice Albedo: Case of Water on Bare Ice
	334	C + -----------------------------------------------------
	335
	336	isn = max(min(isnoSV(ikl) ,ispiSV(ikl)),0)
	337
	338	albWIc = aI1dSV-(aI1dSV-aI2dSV)
	339	. * exp(-(rusnSV(ikl) !
	340	. * (1. -SWS_SV(ikl) ! 0 <=> freezing
	341	. * (1 -min(1,iabs(isn-isnoSV(ikl))))) ! 1 <=> isn=isnoSV
	342	. / ru_dSV)**0.50) !
	343	c albWIc = max(aI1dSV,min(aI2dSV,albWIc+slopSV(ikl)*
	344	c . min(5.,max(1.,dx/10000.))))
	345
	346	SignRo = sign(unun,ro_Ice-5.-ro__SV(ikl,isn)) ! RoSN<920kg/m3
	347	SnowOK = max(zero,SignRo)
	348
	349	albWIc = (1. - SnowOK) * albWIc + SnowOK
	350	. * (aI2dSV + (aI3dSV -aI2dSV)
	351	. * (ro__SV(ikl,isn)-ro_Ice)/(roSdSV-ro_Ice))
	352
	353	c + rocdSV < ro < ro_ice \| aI2dSV< al <aI3dSV (fct of density)
	354	c + ro > ro_ice \| aI1dSV< al <aI2dSV (fct of superficial water content
	355
	356
	357	C +--Integrated Snow/Ice Albedo
	358	C + -------------------------------
	359
	360	a_SII1 = albWIc +(albSn1-albWIc) *SnownH
	361	a_SII1 = min(a_SII1 ,albSn1)
	362
	363	a_SII2 = albWIc +(albSn2-albWIc) *SnownH
	364	a_SII2 = min(a_SII2 ,albSn2)
	365
	366	a_SII3 = albWIc +(albSn3-albWIc) *SnownH
	367	a_SII3 = min(a_SII3 ,albSn3)
[3900]	368
	369	DO i=1,6
	370	a_SII6(i) = albWIc +(albSn6(i)-albWIc) *SnownH
	371	a_SII6(i) = min(a_SII6(i) ,albSn6(i))
	372	ENDDO
	373
[3792]	374	cc #AG agesno = min(agsnSV(ikl,isn) ,AgeMax)
	375	cc #AG a_SII1 = a_SII1 -0.175*agesno/AgeMax
	376	C +... Impurities: Col de Porte Parameter.
	377
	378
	379
	380	C +--Elsewhere Integrated Snow/Ice Albedo
	381	C + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	382	c #cp ELSE
	383	albSII = So1dSV*a_SII1
	384	. + So2dSV*a_SII2
	385	. + So3dSV*a_SII3
	386	c #cp END IF
	387
	388
	389	C +--Integrated Snow/Ice/Soil Albedo
	390	C + -------------------------------
	391
	392	alb1sv(ikl) = albssv(ikl) +(a_SII1-albssv(ikl))*SIcenH
	393	alb1sv(ikl) = min(alb1sv(ikl) ,a_SII1)
	394
	395	alb2sv(ikl) = albssv(ikl) +(a_SII2-albssv(ikl))*SIcenH
	396	alb2sv(ikl) = min(alb2sv(ikl) ,a_SII2)
	397
	398	alb3sv(ikl) = albssv(ikl) +(a_SII3-albssv(ikl))*SIcenH
	399	alb3sv(ikl) = min(alb3sv(ikl) ,a_SII3)
[3900]	400
[3792]	401	albisv(ikl) = albssv(ikl) +(albSII-albssv(ikl))*SIcenH
	402	albisv(ikl) = min(albisv(ikl) ,albSII)
[3900]	403
	404	DO i=1,6
	405	alb6sv(ikl,i) = albssv(ikl) +(a_SII6(i)-albssv(ikl))*SIcenH
	406	alb6sv(ikl,i) = min(alb6sv(ikl,i) ,a_SII6(i))
	407	ENDDO
	408
[3792]	409
	410	C +--Integrated Snow/Ice/Soil Albedo: Clouds Correction! Greuell & all., 1994
	411	C + --------------------------------------------------! Glob.&t Planet.Change
	412	! (9):91-114
	413	alb1sv(ikl) = alb1sv(ikl) + 0.05 (cld_SV(ikl)-0.5)SIcenH
	414	. + dalbed * (1.-cld_SV(ikl))
	415	alb2sv(ikl) = alb2sv(ikl) + 0.05 (cld_SV(ikl)-0.5)SIcenH
	416	. + dalbed * (1.-cld_SV(ikl))
	417	alb3sv(ikl) = alb3sv(ikl) + 0.05 (cld_SV(ikl)-0.5)SIcenH
	418	. + dalbed * (1.-cld_SV(ikl))
[3900]	419	alb6sv(ikl,:) =alb6sv(ikl,:)+0.05 (cld_SV(ikl)-0.5)SIcenH
	420	. + dalbed * (1.-cld_SV(ikl))
[3792]	421	albisv(ikl) = albisv(ikl) + 0.05 (cld_SV(ikl)-0.5)SIcenH
	422	. + dalbed * (1.-cld_SV(ikl))
	423
	424	C +--Integrated Snow/Ice/Soil Albedo: Minimum snow albedo = aI1dSV
	425	C + -------------------------------------------------------------
	426
	427	albedo_old = albisv(ikl)
	428	albisv(ikl) = max(albisv(ikl),aI1dSV * SIcenH
	429	. + albssv(ikl) *(1.0 - SIcenH))
	430	alb1sv(ikl) = alb1sv(ikl) - 1.0/3.0 ! 33 %
	431	. * (albedo_old-albisv(ikl)) / So1dSV
	432	alb2sv(ikl) = alb2sv(ikl) - 1.0/3.0 ! 33 %
	433	. * (albedo_old-albisv(ikl)) / So2dSV
	434	alb3sv(ikl) = alb3sv(ikl) - 1.0/3.0 ! 33 %
	435	. * (albedo_old-albisv(ikl)) / So3dSV
[3900]	436	alb6sv(ikl,1:3) = alb6sv(ikl,1:3) - 1.0/6.0 ! 16 %
	437	. * (albedo_old-albisv(ikl)) / (So1dSV/3)
	438	alb6sv(ikl,4:6) = alb6sv(ikl,4:6) - 1.0/6.0 ! 16 %
	439	. * (albedo_old-albisv(ikl)) / (So2dSV/3)
	440
	441
	442	C +--Integrated Snow/Ice/Soil Albedo: Maximum albedo = 95%
[3792]	443	C + -----------------------------------------------------
	444
	445	albedo_old = albisv(ikl)
	446	albisv(ikl) = min(albisv(ikl),0.95)
	447	alb1sv(ikl) = alb1sv(ikl) - 1.0/3.0 ! 33 %
	448	. * (albedo_old-albisv(ikl)) / So1dSV
	449	alb2sv(ikl) = alb2sv(ikl) - 1.0/3.0 ! 33 %
	450	. * (albedo_old-albisv(ikl)) / So2dSV
	451	alb3sv(ikl) = alb3sv(ikl) - 1.0/3.0 ! 33 %
	452	. * (albedo_old-albisv(ikl)) / So3dSV
[3900]	453	alb6sv(ikl,1:3) = alb6sv(ikl,1:3) - 1.0/6.0 ! 16 %
	454	. * (albedo_old-albisv(ikl)) / (So1dSV/3)
	455	alb6sv(ikl,4:6) = alb6sv(ikl,4:6) - 1.0/6.0 ! 16 %
	456	. * (albedo_old-albisv(ikl)) / (So2dSV/3)
[3792]	457
	458
	459	!Sea Ice/snow permanent-interractive prescription from Nemo
	460	!AO_CK 20/02/2020
	461
	462	! No check if coupling update since MAR and NEMO albedo are too different
	463	! and since MAR albedo is computed on properties that are not in NEMO
	464	! prescription for each time step with NEMO values
	465
	466	c #AO if (LSmask(ikl) .eq. 0 .and. coupling_ao .eq. .true.) then
	467	c #AO if (AOmask(ikl) .eq. 0) then
	468	c #AO albisv(ikl) = (1.-AOmask(ikl))* albAOsisv(ikl)
	469	c #AO. +(AOmask(ikl)*albisv(ikl))
	470	c #AO alb1sv(ikl) = (1.-AOmask(ikl))* albAOsisv(ikl)
	471	c #AO. +(AOmask(ikl)*alb1sv(ikl))
	472	c #AO alb2sv(ikl) = (1.-AOmask(ikl))* albAOsisv(ikl)
	473	c #AO. +(AOmask(ikl)*alb2sv(ikl))
	474	c #AO alb3sv(ikl) = (1.-AOmask(ikl))* albAOsisv(ikl)
	475	c #AO. +(AOmask(ikl)*alb3sv(ikl))
	476	c #AO endif
	477	c #AO endif
	478
	479
	480	alb1sv(ikl) = min(max(zero,alb1sv(ikl)),albmax)
	481	alb2sv(ikl) = min(max(zero,alb2sv(ikl)),albmax)
	482	alb3sv(ikl) = min(max(zero,alb3sv(ikl)),albmax)
[3900]	483
	484	DO i=1,6
	485	alb6sv(ikl,i) = min(max(zero,alb6sv(ikl,i)),albmax)
	486	ENDDO
[3792]	487	END DO
	488
	489
	490	C +--Extinction Coefficient: Exponential Factor
	491	C + ==========================================
	492
	493	DO ikl=1,knonv
	494	sExt_1(ikl) = 1.
	495	sExt_2(ikl) = 1.
	496	sExt_3(ikl) = 1.
	497	sEX_sv(ikl,nsno+1) = 1.
	498
	499	coalb1(ikl) = (1. -alb1sv(ikl))*So1dSV
	500	coalb2(ikl) = (1. -alb2sv(ikl))*So2dSV
	501	coalb3(ikl) = (1. -alb3sv(ikl))*So3dSV
	502	coalbm = coalb1(ikl) +coalb2(ikl) +coalb3(ikl)
	503	coalb1(ikl) = coalb1(ikl) /coalbm
	504	coalb2(ikl) = coalb2(ikl) /coalbm
	505	coalb3(ikl) = coalb3(ikl) /coalbm
	506	END DO
	507
	508	cXF
	509
	510	DO isn=nsno,1,-1
	511	DO ikl=1,knonv
	512	sEX_sv(ikl,isn) = 1.0
	513	!sEX_sv(ikl,isn) = 0.95 ! if MAR is too warm in summer
	514	END DO
	515	END DO
	516
	517	DO ikl=1,knonv
	518	DO isn=max(1,isnoSV(ikl)),1,-1
	519
	520	SignRo = sign(unun, rocdSV - ro__SV(ikl,isn))
	521	SnowOK = max(zero,SignRo) ! Ice Density Threshold
	522
	523	RoFrez = 1.e-3 * ro__SV(ikl,isn) * (1.0-eta_SV(ikl,isn))
	524
	525	OpSqrt = sqrt(max(epsi,SnOpSV(ikl,isn)))
	526	exarg1 = SnowOK 1.e2 max(sbeta1*RoFrez/OpSqrt,sbeta2)
	527	. +(1.0-SnowOK) *sbeta5
	528	exarg2 = SnowOK 1.e2 max(sbeta3*RoFrez/OpSqrt,sbeta4)
	529	. +(1.0-SnowOK) *sbeta5
	530	exarg3 = SnowOK 1.e2 sbeta5
	531	. +(1.0-SnowOK) *sbeta5
	532
	533
	534	C +--Integrated Extinction of Solar Irradiance (Normalized Value)
	535	C + ============================================================
	536
	537	sExt_1(ikl) = sExt_1(ikl)
	538	. * exp(min(0.0,-exarg1 *dzsnSV(ikl,isn)))
	539	sign_0 = sign(unun,eps9 -sExt_1(ikl))
	540	sExt_0 = max(zero,sign_0)*sExt_1(ikl)
	541	sExt_1(ikl) = sExt_1(ikl) -sExt_0
	542
	543	sExt_2(ikl) = sExt_2(ikl)
	544	. * exp(min(0.0,-exarg2 *dzsnSV(ikl,isn)))
	545	sign_0 = sign(unun,eps9 -sExt_2(ikl))
	546	sExt_0 = max(zero,sign_0)*sExt_2(ikl)
	547	sExt_2(ikl) = sExt_2(ikl) -sExt_0
	548
	549	sExt_3(ikl) = sExt_3(ikl)
	550	. * exp(min(0.0,-exarg3 *dzsnSV(ikl,isn)))
	551	sign_0 = sign(unun,eps9 -sExt_3(ikl))
	552	sExt_0 = max(zero,sign_0)*sExt_3(ikl)
	553	sExt_3(ikl) = sExt_3(ikl) -sExt_0
	554
	555	sEX_sv(ikl,isn) = coalb1(ikl) * sExt_1(ikl)
	556	. + coalb2(ikl) * sExt_2(ikl)
	557	. + coalb3(ikl) * sExt_3(ikl)
	558	END DO
	559	END DO
	560
	561	DO isn=0,-nsol,-1
	562	DO ikl=1,knonv
	563	sEX_sv(ikl,isn) = 0.0
	564	END DO
	565	END DO
	566
	567
	568	return
[3900]	569
	570
[3792]	571	end
[3900]	572
	573
	574	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	575	SUBROUTINE albedo_kokhanovsky(lambdamin,lambdamax,
	576	. cossza,dopt,albint)
	577	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	578	! Authors: Hajar El Habchi El Fenniri, Etienne Vignon, Cecile Agosta
	579	! Ghislain Picard
	580	! Routine that calculates the integrated snow spectral albedo between
	581	! lambdamin and lambdamax following Kokhanisky and Zege 2004,
	582	! Scattering optics of snow, Applied Optics, Vol 43, No7
	583	! Code inspired from the snowoptics package of Ghislain Picard:
	584	! https://github.com/ghislainp/snowoptics
	585	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	586
	587
	588	USE VARphy
	589
	590	IMPLICIT NONE
	591
	592	! Inputs
	593	!--------
	594	REAL lambdamin ! minimum wavelength for integration [m]
	595	REAL lambdamax ! maximum wavelength for integration [m]
	596	REAL cossza ! solar zenith angle cosinus
	597	REAL dopt ! optical diameter [m]
	598
	599	!Outputs
	600	!-------
	601	REAL albint
	602
	603	! Local Variables
	604	!-----------------
	605
	606	REAL ropt,cosalb,norm,Pas
	607	REAL SSA,alpha,gamm,R,cos30,alb30
	608	INTEGER i
	609
	610
	611	REAL B_amp ! amplification factor
	612	PARAMETER (B_amp=1.6)
	613
	614	REAL g_asy ! asymetry factor
	615	PARAMETER (g_asy=0.845)
	616
	617	INTEGER nlambda ! length of wavelength vector
	618	PARAMETER(nlambda=200)
	619
	620	REAL lmin
	621	PARAMETER(lmin=185.0E-9)
	622
	623	REAL lmax
	624	PARAMETER(lmax=4000.0E-9)
	625
	626	REAL albmax
	627	PARAMETER(albmax=0.99)
	628
	629	REAL wavelengths(nlambda)
	630	REAL ni(nlambda)
	631
	632	DATA wavelengths / 1.85000000e-07, 2.04170854e-07,
	633	. 2.23341709e-07, 2.42512563e-07,
	634	. 2.61683417e-07, 2.80854271e-07, 3.00025126e-07, 3.19195980e-07,
	635	. 3.38366834e-07, 3.57537688e-07, 3.76708543e-07, 3.95879397e-07,
	636	. 4.15050251e-07, 4.34221106e-07, 4.53391960e-07, 4.72562814e-07,
	637	. 4.91733668e-07, 5.10904523e-07, 5.30075377e-07, 5.49246231e-07,
	638	. 5.68417085e-07, 5.87587940e-07, 6.06758794e-07, 6.25929648e-07,
	639	. 6.45100503e-07, 6.64271357e-07, 6.83442211e-07, 7.02613065e-07,
	640	. 7.21783920e-07, 7.40954774e-07, 7.60125628e-07, 7.79296482e-07,
	641	. 7.98467337e-07, 8.17638191e-07, 8.36809045e-07, 8.55979899e-07,
	642	. 8.75150754e-07, 8.94321608e-07, 9.13492462e-07, 9.32663317e-07,
	643	. 9.51834171e-07, 9.71005025e-07, 9.90175879e-07, 1.00934673e-06,
	644	. 1.02851759e-06, 1.04768844e-06, 1.06685930e-06, 1.08603015e-06,
	645	. 1.10520101e-06, 1.12437186e-06, 1.14354271e-06, 1.16271357e-06,
	646	. 1.18188442e-06, 1.20105528e-06, 1.22022613e-06, 1.23939698e-06,
	647	. 1.25856784e-06, 1.27773869e-06, 1.29690955e-06, 1.31608040e-06,
	648	. 1.33525126e-06, 1.35442211e-06, 1.37359296e-06, 1.39276382e-06,
	649	. 1.41193467e-06, 1.43110553e-06, 1.45027638e-06, 1.46944724e-06,
	650	. 1.48861809e-06, 1.50778894e-06, 1.52695980e-06, 1.54613065e-06,
	651	. 1.56530151e-06, 1.58447236e-06, 1.60364322e-06, 1.62281407e-06,
	652	. 1.64198492e-06, 1.66115578e-06, 1.68032663e-06, 1.69949749e-06,
	653	. 1.71866834e-06, 1.73783920e-06, 1.75701005e-06, 1.77618090e-06,
	654	. 1.79535176e-06, 1.81452261e-06, 1.83369347e-06, 1.85286432e-06,
	655	. 1.87203518e-06, 1.89120603e-06, 1.91037688e-06, 1.92954774e-06,
	656	. 1.94871859e-06, 1.96788945e-06, 1.98706030e-06, 2.00623116e-06,
	657	. 2.02540201e-06, 2.04457286e-06, 2.06374372e-06, 2.08291457e-06,
	658	. 2.10208543e-06, 2.12125628e-06, 2.14042714e-06, 2.15959799e-06,
	659	. 2.17876884e-06, 2.19793970e-06, 2.21711055e-06, 2.23628141e-06,
	660	. 2.25545226e-06, 2.27462312e-06, 2.29379397e-06, 2.31296482e-06,
	661	. 2.33213568e-06, 2.35130653e-06, 2.37047739e-06, 2.38964824e-06,
	662	. 2.40881910e-06, 2.42798995e-06, 2.44716080e-06, 2.46633166e-06,
	663	. 2.48550251e-06, 2.50467337e-06, 2.52384422e-06, 2.54301508e-06,
	664	. 2.56218593e-06, 2.58135678e-06, 2.60052764e-06, 2.61969849e-06,
	665	. 2.63886935e-06, 2.65804020e-06, 2.67721106e-06, 2.69638191e-06,
	666	. 2.71555276e-06, 2.73472362e-06, 2.75389447e-06, 2.77306533e-06,
	667	. 2.79223618e-06, 2.81140704e-06, 2.83057789e-06, 2.84974874e-06,
	668	. 2.86891960e-06, 2.88809045e-06, 2.90726131e-06, 2.92643216e-06,
	669	. 2.94560302e-06, 2.96477387e-06, 2.98394472e-06, 3.00311558e-06,
	670	. 3.02228643e-06, 3.04145729e-06, 3.06062814e-06, 3.07979899e-06,
	671	. 3.09896985e-06, 3.11814070e-06, 3.13731156e-06, 3.15648241e-06,
	672	. 3.17565327e-06, 3.19482412e-06, 3.21399497e-06, 3.23316583e-06,
	673	. 3.25233668e-06, 3.27150754e-06, 3.29067839e-06, 3.30984925e-06,
	674	. 3.32902010e-06, 3.34819095e-06, 3.36736181e-06, 3.38653266e-06,
	675	. 3.40570352e-06, 3.42487437e-06, 3.44404523e-06, 3.46321608e-06,
	676	. 3.48238693e-06, 3.50155779e-06, 3.52072864e-06, 3.53989950e-06,
	677	. 3.55907035e-06, 3.57824121e-06, 3.59741206e-06, 3.61658291e-06,
	678	. 3.63575377e-06, 3.65492462e-06, 3.67409548e-06, 3.69326633e-06,
	679	. 3.71243719e-06, 3.73160804e-06, 3.75077889e-06, 3.76994975e-06,
	680	. 3.78912060e-06, 3.80829146e-06, 3.82746231e-06, 3.84663317e-06,
	681	. 3.86580402e-06, 3.88497487e-06, 3.90414573e-06, 3.92331658e-06,
	682	. 3.94248744e-06, 3.96165829e-06, 3.98082915e-06, 4.00000000e-06/
	683
	684
	685	DATA ni /7.74508407e-10, 7.74508407e-10,
	686	. 7.74508407e-10, 7.74508407e-10,
	687	. 7.74508407e-10, 7.74508407e-10, 7.74508407e-10, 7.74508407e-10,
	688	. 6.98381122e-10, 6.23170274e-10, 5.97655992e-10, 5.84106004e-10,
	689	. 5.44327597e-10, 5.71923510e-10, 6.59723827e-10, 8.05183870e-10,
	690	. 1.03110161e-09, 1.36680386e-09, 1.85161253e-09, 2.56487751e-09,
	691	. 3.56462855e-09, 4.89450926e-09, 6.49252022e-09, 9.62029335e-09,
	692	. 1.32335015e-08, 1.75502184e-08, 2.19240625e-08, 3.03304156e-08,
	693	. 4.07715972e-08, 5.00414911e-08, 7.09722331e-08, 1.00773751e-07,
	694	. 1.31427409e-07, 1.42289041e-07, 1.49066787e-07, 2.01558515e-07,
	695	. 2.99106105e-07, 4.03902086e-07, 4.54292169e-07, 5.21906983e-07,
	696	. 6.27643362e-07, 9.43955678e-07, 1.33464494e-06, 1.97278315e-06,
	697	. 2.31801329e-06, 2.20584676e-06, 1.85568138e-06, 1.73395193e-06,
	698	. 1.73101406e-06, 1.91333936e-06, 2.26413019e-06, 3.23959718e-06,
	699	. 4.94316963e-06, 6.89378896e-06, 1.02237444e-05, 1.21439656e-05,
	700	. 1.31567585e-05, 1.33448288e-05, 1.32000000e-05, 1.31608040e-05,
	701	. 1.33048369e-05, 1.40321464e-05, 1.51526244e-05, 1.80342858e-05,
	702	. 3.82875736e-05, 1.07325259e-04, 2.11961637e-04, 3.82008054e-04,
	703	. 5.30897470e-04, 5.29244735e-04, 4.90876605e-04, 4.33905427e-04,
	704	. 3.77795349e-04, 3.17633099e-04, 2.81078564e-04, 2.57579485e-04,
	705	. 2.42203100e-04, 2.23789060e-04, 2.04306870e-04, 1.87909255e-04,
	706	. 1.73117146e-04, 1.61533186e-04, 1.53420328e-04, 1.47578033e-04,
	707	. 1.42334776e-04, 1.35691466e-04, 1.30495414e-04, 1.36065123e-04,
	708	. 1.70928821e-04, 2.66389730e-04, 4.80957955e-04, 8.25041961e-04,
	709	. 1.21654792e-03, 1.50232875e-03, 1.62316078e-03, 1.61649750e-03,
	710	. 1.53736801e-03, 1.42343711e-03, 1.24459117e-03, 1.02388611e-03,
	711	. 7.89112523e-04, 5.97204264e-04, 4.57152413e-04, 3.62341259e-04,
	712	. 2.99128332e-04, 2.57035569e-04, 2.26992203e-04, 2.07110355e-04,
	713	. 2.05835688e-04, 2.25108810e-04, 2.64262893e-04, 3.23594011e-04,
	714	. 3.93061117e-04, 4.62789970e-04, 5.19664416e-04, 5.59739628e-04,
	715	. 5.93476084e-04, 6.22797885e-04, 6.57484833e-04, 6.92849600e-04,
	716	. 7.26584901e-04, 7.56604648e-04, 7.68009488e-04, 7.65579073e-04,
	717	. 7.50526164e-04, 7.39809972e-04, 7.55622847e-04, 8.05099514e-04,
	718	. 9.67279246e-04, 1.16281559e-03, 1.42570247e-03, 2.04986585e-03,
	719	. 2.93971170e-03, 4.49827711e-03, 7.32537532e-03, 1.18889734e-02,
	720	. 1.85851805e-02, 2.86242532e-02, 4.34131035e-02, 6.37828307e-02,
	721	. 9.24145850e-02, 1.35547945e-01, 1.94143245e-01, 2.54542814e-01,
	722	. 3.02282024e-01, 3.42214181e-01, 3.85475065e-01, 4.38000000e-01,
	723	. 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
	724	. 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
	725	. 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
	726	. 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
	727	. 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
	728	. 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
	729	. 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
	730	. 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
	731	. 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
	732	. 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
	733	. 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
	734	. 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01,
	735	. 4.38000000e-01, 4.38000000e-01, 4.38000000e-01, 4.38000000e-01/
	736
	737
	738	Pas =(lmax-lmin)/nlambda
	739	ropt = dopt/2.0
	740	SSA = 3.0/(rhoIce*ropt)
	741	cos30 = cos(30.0/180.0*pi)
	742
	743
	744	albint=0.
	745	norm=0.
	746
	747	DO i = 1,nlambda
	748	gamm = 4.0 * pi * ni(i) / wavelengths(i)
	749	cosalb = 2.0 / (SSA * rhoice) * B_amp * gamm
	750	alpha = 16. / 3 * cosalb / (1.0 - g_asy)
	751	R = exp(-(alpha*0.5) 3.0 / 7.0 * (1.0 + 2.0 * cossza))
	752	alb30 = exp(-(alpha*0.5) 3.0 / 7.0 * (1.0 + 20 * cos30))
	753
	754	IF ((wavelengths(i).GE.lambdamin).AND.
	755	. (wavelengths(i).LT.lambdamax)) THEN
	756	albint=albint+R*Pas ! rectangle integration -> can be improved with trapezintegration
	757	norm=norm+Pas
	758	ENDIF
	759
	760	END DO
	761
	762	albint=max(0.,min(albint/max(norm,1E-10),albmax))
	763
	764	END
	765

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

Download in other formats:

Original Format