Context Navigation

sisvat_sop.F @ 3946

Last change on this file since 3946 was 3900, checked in by evignon, 4 years ago
Commission de la nouvelle interface LMDZ-SISVAT la prochaine commission consistera a supprimer l'ancien repertoire sisvat et a faire un peu de nettoyage.
File size: 32.0 KB

Line
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 \|
30	C \| DOPsnSV : Snow Optical diameter [m] \|
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
69	USE surface_data, only: iflag_albcalc,correc_alb
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
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
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
139	C +-- wavelength limits [m] for each broad band
140
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/
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))
184
185	C + --For outputs (Etienne)
186	C + ------------------------
187	DOPsnSV(ikl,isn)=SnOpSV(ikl,isn)
188	END DO
189	END DO
190
191
192
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
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	cCA +--Correction of snow albedo for Antarctica/Greenland
212	cCA --------------------------------------------------
213
214
215	albCor = correc_alb
216	c #GL albCor = 1.01
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))
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	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)))
282
283	ENDIF
284
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)
299	albSn6(:) = SnowOKalbSn6(:)+(1.0-SnowOK)max(albSno,minalb)
300
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)
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	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)
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	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))
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	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
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%
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
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	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)
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	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
569
570
571	end
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