Context Navigation

radlwsw.F @ 2187

Last change on this file since 2187 was 1945, checked in by slebonnois, 6 years ago
Correction of a bug in radlwsw Venus
File size: 19.3 KB

Line
1	!
2	! $Header: /home/cvsroot/LMDZ4/libf/phylmd/radlwsw.F,v 1.2 2004/10/27 10:14:46 lmdzadmin Exp $
3	!
4	SUBROUTINE radlwsw(dist, rmu0, fract, zzlev,
5	. paprs, pplay,tsol, t)
6	c
7	c======================================================================
8	c Auteur(s): Z.X. Li (LMD/CNRS) date: 19960719
9	c Objet: interface entre le modele et les rayonnements
10	c Arguments:
11	c dist-----input-R- distance astronomique terre-soleil
12	c rmu0-----input-R- cosinus de l'angle zenithal
13	c fract----input-R- duree d'ensoleillement normalisee
14	c zzlev----input-R- altitude a inter-couche (m)
15	c paprs----input-R- pression a inter-couche (Pa)
16	c pplay----input-R- pression au milieu de couche (Pa)
17	c tsol-----input-R- temperature du sol (en K)
18	c t--------input-R- temperature (K)
19
20	c MODIFS pour multimatrices ksi SPECIFIQUE VENUS
21	c S. Lebonnois 20/12/2006
22	c corrections 13/07/2007
23	c New ksi matrix: possibility of different cloud model fct of lat 05/2014
24
25	c With extension NLTE (G. Gilli, 2014)
26
27	c Ksi matrices latitudinaly interpolated (I. Garate-Lopez, 2016)
28
29	c======================================================================
30	use dimphy
31	USE geometry_mod, ONLY: latitude_deg
32	USE phys_state_var_mod, only: heat,cool,radsol,
33	. topsw,toplw,solsw,sollw,sollwdown,lwnet,swnet
34	use write_field_phy
35	IMPLICIT none
36	#include "YOMCST.h"
37	#include "clesphys.h"
38	#include "comcstVE.h"
39	#include "nlteparams.h"
40
41	!===========
42	! Arguments
43	!===========
44	real rmu0(klon), fract(klon), dist
45
46	REAL zzlev(klon,klev+1)
47	real paprs(klon,klev+1), pplay(klon,klev)
48	real tsol(klon)
49	real t(klon,klev)
50
51	!===========
52	! Local
53	!===========
54	INTEGER k, kk, i, j, band
55
56	REAL PPB(klev+1)
57	REAL PPA(klev)
58
59	REAL zfract, zrmu0,latdeg
60
61	REAL zheat(klev), zcool(klev)
62	real temp(klev),znivs(klev+1)
63	REAL ZFSNET(klev+1),ZFLNET(klev+1)
64	REAL ztopsw, ztoplw
65	REAL zsolsw, zsollw
66	cIM BEG
67	REAL zsollwdown
68	cIM END
69	real,save,allocatable :: ksive(:,:,:,:) ! ksi matrixes in Vincent's file
70
71	real psi(0:klev+1,0:klev+1)
72	real deltapsi(0:klev+1,0:klev+1)
73	real pt0(0:klev+1)
74	real bplck(0:klev+1,nnuve) ! Planck luminances in table layers
75	real y(0:klev,nnuve) ! temporary variable for Planck
76	real zdblay(0:klev+1,nnuve) ! temperature gradient of planck function
77	integer mat0,lat,ips,isza,ips0,isza0
78	real factp,factz,ksi
79	c ------- for lat-interp ----------------
80	integer mat0A, mat0B, latA, latB, kasua
81	integer ipsA, ipsB, iszaA, iszaB, ips0A, ips0B, isza0A, isza0B
82	real lat_deg, latA_deg, latB_deg
83	real factlat, k1, k2, k3, k4
84	c --------------------------------------
85	logical firstcall
86	data firstcall/.true./
87	save firstcall
88
89	cERROR ! For checking if the file it's being read
90	c-------------------------------------------
91	c Initialisations
92	c-----------------
93
94	if (firstcall) then
95
96	c ---------- ksive --------------
97	allocate(ksive(0:klev+1,0:klev+1,nnuve,nbmat))
98	call load_ksi(ksive)
99
100	endif ! firstcall
101	c-------------------------------------------
102
103	DO k = 1, klev
104	DO i = 1, klon
105	heat(i,k)=0.
106	cool(i,k)=0.
107	ENDDO
108	ENDDO
109
110
111	c+++++++ BOUCLE SUR LA GRILLE +++++++++++++++++++++++++
112	DO j = 1, klon
113
114	c======================================================================
115	c Initialisations
116	c ---------------
117
118	DO k = 1, klev
119	zheat(k) = 0.0
120	zcool(k) = 0.0
121	ENDDO
122	c zheat(1:klev)=0.0 !Explicit loop (no change in performance)
123	c zcool(1:klev)=0.0
124
125	DO k = 1, klev+1
126	ZFLNET(k) = 0.0
127	ZFSNET(k) = 0.0
128	ENDDO
129	c ZFLNET(1:klev+1)=0.0
130	c ZFSNET(1:klev+1)=0.0
131
132	ztopsw = 0.0
133	ztoplw = 0.0
134	zsolsw = 0.0
135	zsollw = 0.0
136	zsollwdown = 0.0
137
138	zfract = fract(j)
139	zrmu0 = rmu0(j)
140
141	DO k = 1, klev+1
142	PPB(k) = paprs(j,k)/1.e5
143	ENDDO
144	DO k = 1,klev
145	PPA(k) = pplay(j,k)/1.e5
146	ENDDO
147
148	pt0(0) = tsol(j)
149	DO k = 1, klev
150	pt0(k) = t(j,k)
151	ENDDO
152	pt0(klev+1) = 0.
153
154	DO k = 0,klev+1
155	DO i = 0,klev+1
156	psi(i,k) = 0. ! positif quand nrj de i->k
157	deltapsi(i,k) = 0.
158	ENDDO
159	ENDDO
160
161	c======================================================================
162	c Getting psi and deltapsi
163	c ------------------------
164
165	c Planck function
166	c ---------------
167	do band=1,nnuve
168	do k=0,klev
169	c B(T,l) = al/(exp(bl/T)-1)
170	y(k,band) = exp(bl(band)/pt0(k))-1.
171	bplck(k,band) = al(band)/(y(k,band))
172	zdblay(k,band)= al(band)bl(band)exp(bl(band)/pt0(k))/
173	. ((pt0(k)pt0(k))(y(k,band)*y(k,band)))
174	enddo
175	bplck(klev+1,band) = 0.0
176	zdblay(klev+1,band)= 0.0
177	enddo
178
179	c finding the right matrixes
180	c --------------------------
181
182	mat0 = 0
183	mat0A = 0
184	mat0B = 0
185
186	c Latitude
187	c --------
188	lat = 0
189	latA = 0
190	latB = 0
191
192	c write(,) 'nlatve:', nlatve
193
194	lat_deg = abs(latitude_deg(j))
195
196	c if (nlatve.eq.1) then ! clouds are taken as uniform
197	if ((nlatve.eq.1).or.(lat_deg.le.25.)) then
198	lat = 1
199	elseif (lat_deg.le.50.) then
200	lat = 1
201	latA = 1
202	latB = 2
203	latA_deg = 25.0
204	latB_deg = 55.0
205	elseif (lat_deg.le.55.) then
206	lat = 2
207	latA = 1
208	latB = 2
209	latA_deg = 25.0
210	latB_deg = 55.0
211	elseif (lat_deg.le.60.) then
212	lat = 2
213	latA = 2
214	latB = 3
215	latA_deg = 55.0
216	latB_deg = 65.0
217	elseif (lat_deg.le.65.) then
218	lat = 3
219	latA = 2
220	latB = 3
221	latA_deg = 55.0
222	latB_deg = 65.0
223	elseif (lat_deg.le.70.) then
224	lat = 3
225	latA = 3
226	latB = 4
227	latA_deg = 65.0
228	latB_deg = 75.0
229	elseif (lat_deg.le.75.) then
230	lat = 4
231	latA = 3
232	latB = 4
233	latA_deg = 65.0
234	latB_deg = 75.0
235	elseif (lat_deg.le.80.) then
236	lat = 4
237	latA = 4
238	latB = 5
239	latA_deg = 75.0
240	latB_deg = 85.0
241	elseif (lat_deg.le.85.) then
242	lat = 5
243	latA = 4
244	latB = 5
245	latA_deg = 75.0
246	latB_deg = 85.0
247	else
248	lat = 5
249	endif
250
251	c write(,) 'Lat',lat,'LatA',latA,'LatB',latB
252
253	factlat = 0
254	if (latA.gt.0) then
255	factlat = (lat_deg - latA_deg) / (latB_deg - latA_deg)
256	endif
257
258	c write (,) 'Factor de correccion:', factlat
259
260
261	c Pressure at Surface
262	c -------------------
263
264	ips0=0
265	ips0A=0
266	ips0B=0
267	if (nbpsve(lat).gt.1) then ! Interpolation on ps
268	do ips=1,nbpsve(lat)-1
269	if ( (psurfve(ips,lat).ge.paprs(j,1))
270	. .and.(psurfve(ips+1,lat).lt.paprs(j,1)) ) then
271	ips0 = ips
272	c print*,'ig=',j,' ips0=',ips
273	factp = (paprs(j,1) -psurfve(ips0,lat))
274	. /(psurfve(ips0+1,lat)-psurfve(ips0,lat))
275	exit
276	endif
277	enddo
278	else ! Only one ps, no interpolation
279	ips0=1
280	endif
281
282	if (latA.eq.lat) then
283	ips0A=ips0
284	else
285	if (latA.gt.0) then
286	if (nbpsve(latA).gt.1) then
287	do ipsA=1,nbpsve(latA)-1
288	if ( (psurfve(ipsA,latA).ge.paprs(j,1))
289	. .and.(psurfve(ipsA+1,latA).lt.paprs(j,1)) ) then
290	ips0A = ipsA
291	exit
292	endif
293	enddo
294	else ! Only one ps, no interpolation
295	ips0A=1
296	endif ! nbpsve(latA).gt.1
297	endif ! latA.gt.0 (if latA=0 ips0A is not used, so it doesn't matter)
298	endif ! latA.eq.lat
299
300	if (latB.eq.lat) then
301	ips0B=ips0
302	else
303	if (latB.gt.0) then
304	if (nbpsve(latB).gt.1) then
305	do ipsB=1,nbpsve(latB)-1
306	if ( (psurfve(ipsB,latB).ge.paprs(j,1))
307	. .and.(psurfve(ipsB+1,latB).lt.paprs(j,1)) ) then
308	ips0B = ipsB
309	exit
310	endif
311	enddo
312	else
313	ips0B=1
314	endif ! nbpsve(latB).gt.1
315	endif ! latB.gt.0 (if latB=0 ips0B is not used, so it doesn't matter)
316	endif ! latB.eq.lat
317
318
319	c Solar Zenith Angle
320	c ------------------
321
322	isza0=0
323	isza0A=0
324	isza0B=0
325	if (nbszave(lat).gt.1) then
326	do isza=1,nbszave(lat)-1
327	if ( (szave(isza,lat).ge.zrmu0)
328	. .and.(szave(isza+1,lat).lt.zrmu0) ) then
329	isza0 = isza
330	c print*,'ig=',j,' isza0=',isza
331	factz = (zrmu0 -szave(isza0,lat))
332	. /(szave(isza0+1,lat)-szave(isza0,lat))
333	exit
334	endif
335	enddo
336	else ! Only one sza, no interpolation
337	isza0=-99
338	endif
339
340
341	if (latA.eq.lat) then
342	isza0A=isza0
343	else
344	if (latA.gt.0) then
345	if (nbszave(latA).gt.1) then
346	do iszaA=1,nbszave(latA)-1
347	if ( (szave(iszaA,latA).ge.zrmu0)
348	. .and.(szave(iszaA+1,latA).lt.zrmu0) ) then
349	isza0A = iszaA
350	exit
351	endif
352	enddo
353	else ! Only one sza, no interpolation
354	isza0A=-99
355	endif ! nbszave(latA).gt.1
356	endif ! latA.gt.0 (if latA=0 isza0A is not used, so it doesn't matter)
357	endif ! latA.eq.lat
358
359	if (latB.eq.lat) then
360	isza0B=isza0
361	else
362	if (latB.gt.0) then
363	if (nbszave(latB).gt.1) then
364	! init to avoid outside values (near midnight so similar compo...)
365	isza0B = nbszave(latB)
366	do iszaB=1,nbszave(latB)-1
367	if ( (szave(iszaB,latB).ge.zrmu0)
368	. .and.(szave(iszaB+1,latB).lt.zrmu0) ) then
369	isza0B = iszaB
370	exit
371	endif
372	enddo
373	else ! Only one sza, no interpolation
374	isza0B=-99
375	endif ! nbszave(latB).gt.1
376	endif ! latB.gt.0 (if latB=0 isza0B is not used, so it doesn't matter)
377	endif ! latB.eq.lat
378
379	c write(,) 'nbszave', nbszave(lat),'nbpsve(lat)',nbpsve(lat)
380
381
382	c -------- Probleme aux bords
383	c surf press lower than the lowest surf pres in matrices
384	if ((ips0.eq.0).and.(psurfve(nbpsve(lat),lat).gt.paprs(j,1)))
385	. then
386	ips0 = nbpsve(lat)-1
387	print*,'Extrapolation! ig=',j,' ips0=',ips0
388	factp = (paprs(j,1) -psurfve(ips0,lat))
389	. /(psurfve(ips0+1,lat)-psurfve(ips0,lat))
390	endif
391	c surf press higher than the highest surf pres in matrices
392	if ((ips0.eq.0).and.(psurfve(1,lat).le.paprs(j,1))) then
393	ips0 = 1
394	print*,'Extrapolation! ig=',j,' ips0=',ips0
395	factp = (paprs(j,1) -psurfve(ips0,lat))
396	. /(psurfve(ips0+1,lat)-psurfve(ips0,lat))
397	endif
398
399	c this has to be done for ips0A and ips0B also...
400	if (latA.eq.lat) then
401	ips0A = ips0
402	else
403	if (latA.gt.0) then
404	if ((ips0A.eq.0).and.
405	. (psurfve(nbpsve(latA),latA).gt.paprs(j,1))) then
406	ips0A = nbpsve(latA)-1
407	endif
408	if ((ips0A.eq.0).and.(psurfve(1,latA).le.paprs(j,1))) then
409	ips0A = 1
410	endif
411	endif
412	endif
413	if (latB.eq.lat) then
414	ips0B = ips0
415	else
416	if (latB.gt.0) then
417	if ((ips0B.eq.0).and.
418	. (psurfve(nbpsve(latB),latB).gt.paprs(j,1))) then
419	ips0B = nbpsve(latB)-1
420	endif
421	if ((ips0B.eq.0).and.(psurfve(1,latB).le.paprs(j,1))) then
422	ips0B = 1
423	endif
424	endif
425	endif
426
427	c ---------
428
429	if ((ips0.eq.0).or.(isza0.eq.0)) then
430	write(,) 'Finding the right matrix in radlwsw'
431	print*,'Interpolation problem, grid point ig=',j
432	print*,'psurf = ',paprs(j,1),' mu0 = ',zrmu0
433	stop
434	endif
435
436	if (isza0.eq.-99) then
437	mat0 = indexve(lat) +ips0
438	if (latA.gt.0) then
439	mat0A = indexve(latA)+ips0A
440	mat0B = indexve(latB)+ips0B
441	endif
442	else
443	mat0 = indexve(lat) +(isza0 -1)*nbpsve(lat) +ips0
444	if (latA.gt.0) then
445	mat0A = indexve(latA)+(isza0A-1)*nbpsve(latA)+ips0A
446	mat0B = indexve(latB)+(isza0B-1)*nbpsve(latB)+ips0B
447	endif
448	endif
449
450	c write(,) 'Second revision> Lat',lat,'LatA',latA,'LatB',latB
451
452	c interpolation of ksi and computation of psi,deltapsi
453	c ----------------------------------------------------
454
455	if (isza0.eq.-99) then
456	if (latA.gt.0) then ! Not being in the two extremal bins
457
458	do band=1,nnuve
459	do k=0,klev+1
460	do i=k+1,klev+1
461	k1 = ksive(i,k,band,mat0A)*(1-factlat)
462	. + ksive(i,k,band,mat0B)*factlat
463	k2 = ksive(i,k,band,mat0A+1)*(1-factlat)
464	. + ksive(i,k,band,mat0B+1)*factlat
465	ksi = k1(1-factp) + k2factp
466	psi(i,k) = psi(i,k) +
467	. RPIksi(bplck(i,band)-bplck(k,band))
468	c ONLY NEEDED IF IMPLICIT CHOSEN IN LW_VENUS_VE (not the case right now)
469	c deltapsi(i,k) = deltapsi(i,k) + RPIksizdblay(i,band)
470	c deltapsi(k,i) = deltapsi(k,i) + RPIksizdblay(k,band)
471
472	kasua=1
473	enddo
474	enddo
475	enddo
476	do k=0,klev+1
477	do i=k+1,klev+1
478	psi(k,i) = -psi(i,k)
479	enddo
480	enddo
481
482	else ! latA=0 --> extremal bins
483
484	do band=1,nnuve
485	do k=0,klev+1
486	do i=k+1,klev+1
487	ksi = ksive(i,k,band,mat0)*(1-factp)
488	. + ksive(i,k,band,mat0+1)*factp
489	psi(i,k) = psi(i,k) +
490	. RPIksi(bplck(i,band)-bplck(k,band))
491	c ONLY NEEDED IF IMPLICIT CHOSEN IN LW_VENUS_VE (not the case right now)
492	c deltapsi(i,k) = deltapsi(i,k) + RPIksizdblay(i,band)
493	c deltapsi(k,i) = deltapsi(k,i) + RPIksizdblay(k,band)
494
495	kasua=2
496	enddo
497	enddo
498	enddo
499	do k=0,klev+1
500	do i=k+1,klev+1
501	psi(k,i) = -psi(i,k)
502	enddo
503	enddo
504
505	endif ! latA.gt.0
506
507	else ! isza0=!-99
508
509	if (latA.gt.0) then ! Not being in the two extremal bins
510
511	do band=1,nnuve
512	do k=0,klev+1
513	do i=k+1,klev+1
514	k1 = ksive(i,k,band,mat0A)*(1-factlat)
515	. + ksive(i,k,band,mat0B)*factlat
516	k2 = ksive(i,k,band,mat0A+1)*(1-factlat)
517	. + ksive(i,k,band,mat0B+1)*factlat
518	k3 = ksive(i,k,band,mat0A+nbpsve(latA))*(1-factlat)
519	. + ksive(i,k,band,mat0B+nbpsve(latB))*factlat
520	k4 = ksive(i,k,band,mat0A+nbpsve(latA)+1)*(1-factlat)
521	. + ksive(i,k,band,mat0B+nbpsve(latB)+1)*factlat
522	ksi = ( k1(1-factp) + k2factp )*(1-factz)
523	. + ( k3(1-factp) + k4factp )*factz
524	psi(i,k) = psi(i,k) +
525	. RPIksi(bplck(i,band)-bplck(k,band))
526	c ONLY NEEDED IF IMPLICIT CHOSEN IN LW_VENUS_VE (not the case right now)
527	c deltapsi(i,k) = deltapsi(i,k) + RPIksizdblay(i,band)
528	c deltapsi(k,i) = deltapsi(k,i) + RPIksizdblay(k,band)
529
530	kasua=3
531	enddo
532	enddo
533	enddo
534	do k=0,klev+1
535	do i=k+1,klev+1
536	psi(k,i) = -psi(i,k)
537	enddo
538	enddo
539
540	else ! latA=0 --> extremal bins
541
542	do band=1,nnuve
543	do k=0,klev+1
544	do i=k+1,klev+1
545	ksi = ksive(i,k,band,mat0)(1-factp)(1-factz)
546	. + ksive(i,k,band,mat0+1)factp (1-factz)
547	. + ksive(i,k,band,mat0+nbpsve(lat))(1-factp)factz
548	. + ksive(i,k,band,mat0+nbpsve(lat)+1)factp factz
549	psi(i,k) = psi(i,k) +
550	. RPIksi(bplck(i,band)-bplck(k,band))
551	c ONLY NEEDED IF IMPLICIT CHOSEN IN LW_VENUS_VE (not the case right now)
552	c deltapsi(i,k) = deltapsi(i,k) + RPIksizdblay(i,band)
553	c deltapsi(k,i) = deltapsi(k,i) + RPIksizdblay(k,band)
554
555	kasua=4
556	enddo
557	enddo
558	enddo
559	do k=0,klev+1
560	do i=k+1,klev+1
561	psi(k,i) = -psi(i,k)
562	enddo
563	enddo
564
565	endif ! latA.gt.0
566	endif ! isza0.eq.-99
567
568	c write(,) 'Kasua:', kasua
569
570	c======================================================================
571	c LW call
572	c---------
573	temp(1:klev)=t(j,1:klev)
574	CALL LW_venus_ve(
575	. PPB,temp,psi,deltapsi,
576	. zcool,
577	. ztoplw,zsollw,
578	. zsollwdown,ZFLNET)
579	c---------
580	c SW call
581	c---------
582	znivs=zzlev(j,:)
583	latdeg=abs(latitude_deg(j))
584
585	c CALL SW_venus_ve_1Dglobave(zrmu0, zfract, ! pour moy globale
586	c CALL SW_venus_ve(zrmu0, zfract,
587	c S PPB,temp,znivs,
588	c S zheat,
589	c S ztopsw,zsolsw,ZFSNET)
590
591	c CALL SW_venus_cl_1Dglobave(zrmu0,zfract, ! pour moy globale
592	c CALL SW_venus_cl(zrmu0,zfract,
593	c CALL SW_venus_dc_1Dglobave(zrmu0,zfract, ! pour moy globale
594	c CALL SW_venus_dc(zrmu0,zfract,
595	c CALL SW_venus_rh_1Dglobave(zrmu0,zfract, ! pour moy globale
596	c S PPB,temp,
597	CALL SW_venus_rh(zrmu0,zfract,latdeg,
598	S PPA,PPB,temp,
599	S zheat,
600	S ztopsw,zsolsw,ZFSNET)
601	c======================================================================
602	radsol(j) = zsolsw - zsollw ! + vers bas
603	topsw(j) = ztopsw ! + vers bas
604	toplw(j) = ztoplw ! + vers haut
605	solsw(j) = zsolsw ! + vers bas
606	sollw(j) = -zsollw ! + vers bas
607	sollwdown(j) = zsollwdown ! + vers bas
608
609	DO k = 1, klev+1
610	lwnet (j,k) = ZFLNET(k)
611	swnet (j,k) = ZFSNET(k)
612	ENDDO
613
614	c
615	C heat/cool with upper atmosphere
616	C
617	IF(callnlte) THEN
618	DO k = 1,nlaylte
619	heat(j,k) = zheat(k)
620	cool(j,k) = zcool(k)
621	ENDDO
622	c Zero tendencies for any remaining layers between nlaylte and klev
623	if (klev.gt.nlaylte) then
624	do k = nlaylte+1, klev
625	heat(j,k) = 0.
626	cool(j,k) = 0.
627	enddo
628	endif
629	ELSE
630	DO k = 1, klev
631	heat(j,k) = zheat(k)
632	cool(j,k) = zcool(k)
633	ENDDO
634	ENDIF ! callnlte
635	ENDDO ! of DO j = 1, klon
636	c+++++++ FIN BOUCLE SUR LA GRILLE +++++++++++++++++++++++++
637
638	! for tests: write output fields...
639	! call writefield_phy('radlwsw_heat',heat,klev)
640	! call writefield_phy('radlwsw_cool',cool,klev)
641	! call writefield_phy('radlwsw_radsol',radsol,1)
642	! call writefield_phy('radlwsw_topsw',topsw,1)
643	! call writefield_phy('radlwsw_toplw',toplw,1)
644	! call writefield_phy('radlwsw_solsw',solsw,1)
645	! call writefield_phy('radlwsw_sollw',sollw,1)
646	! call writefield_phy('radlwsw_sollwdown',sollwdown,1)
647	! call writefield_phy('radlwsw_swnet',swnet,klev+1)
648	! call writefield_phy('radlwsw_lwnet',lwnet,klev+1)
649
650	c tests
651
652	c j = klon/2
653	c j = 1
654	c print*,'mu0=',rmu0(j)
655	c print*,' net flux vis HEAT(K/Eday)'
656	c do k=1,klev
657	c print,k,ZFSNET(k),heat(j,k)86400.
658	c enddo
659	c print*,' net flux IR COOL(K/Eday)'
660	c do k=1,klev
661	c print,k,ZFLNET(k),cool(j,k)86400.
662	c enddo
663
664	firstcall = .false.
665	RETURN
666	END

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

Download in other formats:

Original Format