Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

sisvat_sno_albedo.f90 @ 5407

Last change on this file since 5407 was 5246, checked in by abarral, 8 weeks ago
Convert fixed-form to free-form sources .F -> .{f,F}90 (WIP: some .F remain, will be handled in subsequent commits)
File size: 30.9 KB

Rev	Line
[5246]	1	subroutine SnOptP(jjtime)
[3792]	2
[5246]	3	! +------------------------------------------------------------------------+
	4	! \| MAR/SISVAT SnOptP 12-08-2019 MAR \|
	5	! \| SubRoutine SnOptP computes the Snow Pack optical Properties \|
	6	! +------------------------------------------------------------------------+
	7	! \| \|
	8	! \| PARAMETERS: knonv: Total Number of columns = \|
	9	! \| ^^^^^^^^^^ = Total Number of continental Grid Boxes \|
	10	! \| X Number of Mosaic Cell per Grid Box \|
	11	! \| \|
	12	! \| INPUT: isnoSV = total Nb of Ice/Snow Layers \|
	13	! \| ^^^^^ ispiSV = 0,...,nsno: Uppermost Superimposed Ice Layer \|
	14	! \| \|
	15	! \| \|
	16	! \| INPUT: G1snSV : Dendricity (<0) or Sphericity (>0) of Snow Layer \|
	17	! \| ^^^^^ G2snSV : Sphericity (>0) or Size of Snow Layer \|
	18	! \| agsnSV : Snow Age [day] \|
	19	! \| ro__SV : Snow/Soil Volumic Mass [kg/m3] \|
	20	! \| eta_SV : Water Content [m3/m3] \|
	21	! \| rusnSV : Surficial Water Thickness [kg/m2] .OR. [mm] \|
	22	! \| SWS_SV : Surficial Water Status \|
	23	! \| dzsnSV : Snow Layer Thickness [m] \|
	24	! \| \|
	25	! \| albssv : Soil Albedo [-] \|
	26	! \| zzsnsv : Snow Pack Thickness [m] \|
	27	! \| \|
	28	! \| OUTPUT: albisv : Snow/Ice/Water/Soil Integrated Albedo [-] \|
	29	! \| ^^^^^^ sEX_sv : Verticaly Integrated Extinction Coefficient \|
	30	! \| DOPsnSV : Snow Optical diameter [m] \|
	31	! \| \|
	32	! \| Internal Variables: \|
	33	! \| ^^^^^^^^^^^^^^^^^^ \|
	34	! \| SnOpSV : Snow Grain optical Size [m] \|
	35	! \| EX1_sv : Integrated Snow Extinction (0.3--0.8micr.m) \|
	36	! \| EX2_sv : Integrated Snow Extinction (0.8--1.5micr.m) \|
	37	! \| EX3_sv : Integrated Snow Extinction (1.5--2.8micr.m) \|
	38	! \| \|
	39	! \| METHODE: Calcul de la taille optique des grains ? partir de \|
	40	! \| ^^^^^^^ -leur type decrit par les deux variables descriptives \|
	41	! \| continues sur la plage -99/+99 passees en appel. \|
	42	! \| -la taille optique (1/10mm) des etoiles, \|
	43	! \| des grains fins et \|
	44	! \| des jeunes faces planes \|
	45	! \| \|
	46	! \| METHOD: Computation of the optical diameter of the grains \|
	47	! \| ^^^^^^ described with the CROCUS formalism G1snSV / G2snSV \|
	48	! \| \|
	49	! \| REFERENCE: Brun et al. 1989, J. Glaciol 35 pp. 333--342 \|
	50	! \| ^^^^^^^^^ Brun et al. 1992, J. Glaciol 38 pp. 13-- 22 \|
	51	! \| Eric Martin Sept.1996 \|
	52	! \| \|
	53	! \| \|
	54	! +------------------------------------------------------------------------+
[3792]	55
	56
	57
	58
[5246]	59	! +--Global Variables
	60	! + ================
[3792]	61
	62
[5246]	63	use VARphy
	64	use VAR_SV
	65	use VARdSV
	66	use VARxSV
	67	use VARySV
	68	use VARtSV
	69	USE surface_data, only: iflag_albcalc,correc_alb
[3900]	70
[5246]	71	IMPLICIT NONE
[3900]	72
	73
[5246]	74	! + -- INPUT
	75	integer :: jjtime
[3900]	76
[5246]	77	! +--Internal Variables
	78	! + ==================
[3792]	79
[5246]	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	! #AG real agesno
[3900]	89
[5246]	90	integer :: isn ,ikl ,isn1, i
	91	real :: sbeta1,sbeta2,sbeta3,sbeta4,sbeta5
	92	real :: AgeMax,AlbMin,HSnoSV,HIceSV,doptmx,SignG1,Sph_OK
	93	real :: dalbed,dalbeS,dalbeW
	94	real :: bsegal,czemax,csegal,csza
	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
	104	real :: l1min,l1max,l2min,l2max,l3min,l3max
	105	real :: l6min(6), l6max(6), albSn6(6), a_SII6(6)
	106	real :: lmintmp,lmaxtmp,albtmp
[3900]	107
[5246]	108	! +--Local DATA
	109	! + ============
[3900]	110
[5246]	111	! +--For the computation of the solar irradiance extinction in snow
	112	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	113	data sbeta1/0.0192/,sbeta2/0.4000/,sbeta3/0.1098/
	114	data sbeta4/1.0000/
	115	data sbeta5/2.00e1/
[3900]	116
[5246]	117	! +--Snow Age Maximum (Taiga, e.g. Col de Porte)
	118	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	119	data AgeMax /60.0/
	120	! +... AgeMax: Snow Age Maximum [day]
[3900]	121
[5246]	122	data AlbMin /0.94/
	123	! +... AlbMin: Albedo Minimum / visible (0.3--0.8 micrometers)
[3900]	124
[5246]	125	data HSnoSV /0.01/
	126	! +... HSnoSV: Snow Thickness through witch
	127	! + Albedo is interpolated to Ice Albedo
	128	data HIceSV /0.10/
	129	! +... HIceSV: Snow/Ice Thickness through witch
	130	! + Albedo is interpolated to Soil Albedo
	131	data doptmx /2.3e-3/
	132	! +... doptmx: Maximum optical Diameter (pi * R**2) [m]
	133	! +
	134	data czeMAX /0.173648178/ ! 80.deg (Segal et al., 1991 JAS)
[3900]	135
[5246]	136	data bsegal /4.00 / !
	137	data albmax /0.99 / ! Albedo max
[3900]	138
[5246]	139	! +-- wavelength limits [m] for each broad band
[3900]	140
[5246]	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/
[3900]	144
[5246]	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/
[3900]	149
	150
[5246]	151	! +--Snow Grain optical Size
	152	! + =======================
[3900]	153
[5246]	154	DO ikl=1,knonv
	155	DO isn=1,max(1,isnoSV(ikl))
[3900]	156
[5246]	157	G2snSV(ikl,isn) = max(epsi,G2snSV(ikl,isn))
	158	! +... Avoid non physical Values
[3900]	159
[5246]	160	SignG1 = sign(unun,G1snSV(ikl,isn))
	161	Sph_OK = max(zero,SignG1)
[3900]	162
[5246]	163	SnOpSV(ikl,isn) = 1.e-4 * &
	164	! +... SI: (from 1/10 mm to m)
[3900]	165
[3792]	166
[5246]	167	! +--Contribution of Non Dendritic Snow
	168	! + ----------------------------------
[3900]	169
[5246]	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)) &
[3900]	173
	174
[5246]	175	! +--Contribution of Dendritic Snow
	176	! + ----------------------------------
[3900]	177
[5246]	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
[5246]	185	! + --For outputs (Etienne)
	186	! + ------------------------
	187	DOPsnSV(ikl,isn)=SnOpSV(ikl,isn)
	188	END DO
	189	END DO
[3792]	190
	191
	192
	193
[5246]	194	! +--Snow/Ice Albedo
	195	! + ===============
	196
	197
	198
	199	! +--Uppermost effective Snow Layer
	200	! + ------------------------------
	201
	202	DO ikl=1,knonv
	203
	204	isn = max(iun,isnoSV(ikl))
	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	!CA +--Correction of snow albedo for Antarctica/Greenland
	212	!CA --------------------------------------------------
	213
	214
	215	albCor = correc_alb
	216	! #GL albCor = 1.01
	217	! #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	! #cz dalbeS = ((bsegal+unun)/(unun+2.0bsegalcsegal)
	249	! #cz. - unun )*0.32
	250	! #cz. / bsegal
	251	! #cz dalbeS = max(dalbeS,zero)
	252	! #cz dalbed = dalbeS * min(1,isnoSV(ikl))
	253
	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
	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
	275	albSn6(1:3)=albSn1
	276	albSn6(4:6)=albSn2
	277
	278	! !snow albedo corection if wetsnow
	279	! #GL albSn1 = albSn1max(0.9,(1.-1.5eta_SV(ikl,isn)))
	280	! #GL albSn2 = albSn2max(0.9,(1.-1.5eta_SV(ikl,isn)))
	281	! #GL albSn3 = albSn3max(0.9,(1.-1.5eta_SV(ikl,isn)))
	282
	283	ENDIF
	284
	285
	286	albSno = So1dSV*albSn1 &
	287	+ So2dSV*albSn2 &
	288	+ So3dSV*albSn3
	289
	290	!XF
	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)
	299	albSn6(:) = SnowOKalbSn6(:)+(1.0-SnowOK)max(albSno,minalb)
	300
	301
	302	! + ro < roSdSV \| min al > aI3dSV
	303	! + roSdSV < ro < rocdSV \| aI2dSV < min al < aI3dSV (fct of density)
	304
	305
	306	! +--Snow/Ice Pack Thickness
	307	! + -----------------------
	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	! + The value of SnownH is set to 1 in case of ice lenses above
	318	! + 1m of dry snow (ro<600kg/m3) for using CROCUS albedo
	319
	320	! ro_ave = 0.
	321	! dz_ave = 0.
	322	! SnowOK = 1.
	323	! do isn = isnoSV(ikl),1,-1
	324	! ro_ave = ro_ave + ro__SV(ikl,isn) * dzsnSV(ikl,isn) * SnowOK
	325	! dz_ave = dz_ave + dzsnSV(ikl,isn) * SnowOK
	326	! SnowOK = max(zero,sign(unun,1.-dz_ave))
	327	! enddo
	328
	329	! ro_ave = ro_ave / max(dz_ave,epsi)
	330	! SnowOK = max(zero,sign(unun,600.-ro_ave))
	331	! SnownH = SnowOK + SnownH * (1. - SnowOK)
	332
	333	! +--Integrated Snow/Ice Albedo: Case of Water on Bare Ice
	334	! + -----------------------------------------------------
	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	! albWIc = max(aI1dSV,min(aI2dSV,albWIc+slopSV(ikl)*
	344	! . 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	! + rocdSV < ro < ro_ice \| aI2dSV< al <aI3dSV (fct of density)
	354	! + ro > ro_ice \| aI1dSV< al <aI2dSV (fct of superficial water content
	355
	356
	357	! +--Integrated Snow/Ice Albedo
	358	! + -------------------------------
	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)
	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
	374	!c #AG agesno = min(agsnSV(ikl,isn) ,AgeMax)
	375	!c #AG a_SII1 = a_SII1 -0.175*agesno/AgeMax
	376	! +... Impurities: Col de Porte Parameter.
	377
	378
	379
	380	! +--Elsewhere Integrated Snow/Ice Albedo
	381	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	382	! #cp ELSE
	383	albSII = So1dSV*a_SII1 &
	384	+ So2dSV*a_SII2 &
	385	+ So3dSV*a_SII3
	386	! #cp END IF
	387
	388
	389	! +--Integrated Snow/Ice/Soil Albedo
	390	! + -------------------------------
	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)
	400
	401	albisv(ikl) = albssv(ikl) +(albSII-albssv(ikl))*SIcenH
	402	albisv(ikl) = min(albisv(ikl) ,albSII)
	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
	409
	410	! +--Integrated Snow/Ice/Soil Albedo: Clouds Correction! Greuell & all., 1994
	411	! + --------------------------------------------------! 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))
	419	alb6sv(ikl,:) =alb6sv(ikl,:)+0.05 (cld_SV(ikl)-0.5)SIcenH &
	420	+ dalbed * (1.-cld_SV(ikl))
	421	albisv(ikl) = albisv(ikl) + 0.05 (cld_SV(ikl)-0.5)SIcenH &
	422	+ dalbed * (1.-cld_SV(ikl))
	423
	424	! +--Integrated Snow/Ice/Soil Albedo: Minimum snow albedo = aI1dSV
	425	! + -------------------------------------------------------------
	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
	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	! +--Integrated Snow/Ice/Soil Albedo: Maximum albedo = 95%
	443	! + -----------------------------------------------------
	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
	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)
	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	! #AO if (LSmask(ikl) .eq. 0 .and. coupling_ao .eq. .true.) then
	467	! #AO if (AOmask(ikl) .eq. 0) then
	468	! #AO albisv(ikl) = (1.-AOmask(ikl))* albAOsisv(ikl)
	469	! #AO. +(AOmask(ikl)*albisv(ikl))
	470	! #AO alb1sv(ikl) = (1.-AOmask(ikl))* albAOsisv(ikl)
	471	! #AO. +(AOmask(ikl)*alb1sv(ikl))
	472	! #AO alb2sv(ikl) = (1.-AOmask(ikl))* albAOsisv(ikl)
	473	! #AO. +(AOmask(ikl)*alb2sv(ikl))
	474	! #AO alb3sv(ikl) = (1.-AOmask(ikl))* albAOsisv(ikl)
	475	! #AO. +(AOmask(ikl)*alb3sv(ikl))
	476	! #AO endif
	477	! #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)
	483
	484	DO i=1,6
	485	alb6sv(ikl,i) = min(max(zero,alb6sv(ikl,i)),albmax)
	486	ENDDO
	487	END DO
	488
	489
	490	! +--Extinction Coefficient: Exponential Factor
	491	! + ==========================================
	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	!XF
	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
[3792]	514	END DO
[5246]	515	END DO
[3900]	516
[5246]	517	DO ikl=1,knonv
	518	DO isn=max(1,isnoSV(ikl)),1,-1
[3900]	519
[5246]	520	SignRo = sign(unun, rocdSV - ro__SV(ikl,isn))
	521	SnowOK = max(zero,SignRo) ! Ice Density Threshold
[3900]	522
[5246]	523	RoFrez = 1.e-3 * ro__SV(ikl,isn) * (1.0-eta_SV(ikl,isn))
[3900]	524
[5246]	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	! +--Integrated Extinction of Solar Irradiance (Normalized Value)
	535	! + ============================================================
	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
	569
	570
	571	end subroutine snoptp
	572
	573
[3900]	574	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[5246]	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	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[3900]	586
	587
[5246]	588	USE VARphy
[3900]	589
[5246]	590	IMPLICIT NONE
[3900]	591
[5246]	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]
[3900]	598
[5246]	599	!Outputs
	600	!-------
	601	REAL :: albint
[3900]	602
[5246]	603	! Local Variables
	604	!-----------------
[3900]	605
[5246]	606	REAL :: ropt,cosalb,norm,Pas
	607	REAL :: SSA,alpha,gamm,R,cos30,alb30
	608	INTEGER :: i
[3900]	609
	610
[5246]	611	REAL :: B_amp ! amplification factor
	612	PARAMETER (B_amp=1.6)
[3900]	613
[5246]	614	REAL :: g_asy ! asymetry factor
	615	PARAMETER (g_asy=0.845)
[3900]	616
[5246]	617	INTEGER :: nlambda ! length of wavelength vector
	618	PARAMETER(nlambda=200)
[3900]	619
[5246]	620	REAL :: lmin
	621	PARAMETER(lmin=185.0E-9)
[3900]	622
[5246]	623	REAL :: lmax
	624	PARAMETER(lmax=4000.0E-9)
[3900]	625
[5246]	626	REAL :: albmax
	627	PARAMETER(albmax=0.99)
[3900]	628
[5246]	629	REAL :: wavelengths(nlambda)
	630	REAL :: ni(nlambda)
[3900]	631
[5246]	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/
[3900]	683
	684
[5246]	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/
[3900]	736
	737
[5246]	738	Pas =(lmax-lmin)/nlambda
	739	ropt = dopt/2.0
	740	SSA = 3.0/(rhoIce*ropt)
	741	cos30 = cos(30.0/180.0*pi)
[3900]	742
	743
[5246]	744	albint=0.
	745	norm=0.
[3900]	746
[5246]	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))
[3900]	753
[5246]	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
[3900]	759
[5246]	760	END DO
	761
	762	albint=max(0.,min(albint/max(norm,1E-10),albmax))
	763
	764	END SUBROUTINE albedo_kokhanovsky
	765

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

Download in other formats:

Original Format