Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sisvat_sno_albedo.f90 @ 5354

Last change on this file since 5354 was 5246, checked in by abarral, 7 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

Line
1	subroutine SnOptP(jjtime)
2
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	! +------------------------------------------------------------------------+
55
56
57
58
59	! +--Global Variables
60	! + ================
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	! + -- INPUT
75	integer :: jjtime
76
77	! +--Internal Variables
78	! + ==================
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	! #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	! +--Local DATA
109	! + ============
110
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/
116
117	! +--Snow Age Maximum (Taiga, e.g. Col de Porte)
118	! + ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
119	data AgeMax /60.0/
120	! +... AgeMax: Snow Age Maximum [day]
121
122	data AlbMin /0.94/
123	! +... AlbMin: Albedo Minimum / visible (0.3--0.8 micrometers)
124
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)
135
136	data bsegal /4.00 / !
137	data albmax /0.99 / ! Albedo max
138
139	! +-- 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	! +--Snow Grain optical Size
152	! + =======================
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	! +... 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	! +... SI: (from 1/10 mm to m)
165
166
167	! +--Contribution of Non Dendritic Snow
168	! + ----------------------------------
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	! +--Contribution of Dendritic Snow
176	! + ----------------------------------
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	! + --For outputs (Etienne)
186	! + ------------------------
187	DOPsnSV(ikl,isn)=SnOpSV(ikl,isn)
188	END DO
189	END DO
190
191
192
193
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
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	! +--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
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 SUBROUTINE albedo_kokhanovsky
765

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

Download in other formats:

Original Format