Context Navigation

xvik.F @ 2397

Last change on this file since 2397 was 2325, checked in by tpierron, 5 years ago

Mars GCM :

Update of the program xvik.F
Useless code removed
Add some comments about xvik outputs
Declaration of physics constants inside the program as variable and not as inputs like before
Temporal shift due to viking landing sol removed : xvik outputs are given in real sol
Add the possibility to choose the format of time dimension for xvik outputs : ls, sol or both
Path to input and output files removed from the code and put as input
List of input files removed from code and put as input

File size: 20.3 KB

Line
1	PROGRAM xvik
2
3	USE filtreg_mod, ONLY: inifilr
4	USE comconst_mod, ONLY: dtvr,g,r,pi
5
6
7	IMPLICIT NONE
8
9
10	c=======================================================================
11	c
12	c Pression au site Viking
13	c
14	c=======================================================================
15
16
17	c-----------------------------------------------------------------------
18	c declarations:
19	c-----------------------------------------------------------------------
20
21
22	include "dimensions.h"
23	include "paramet.h"
24	include "comdissip.h"
25	include "comgeom2.h"
26	include "netcdf.inc"
27
28
29	INTEGER itau,nbpas,nbpasmx
30	PARAMETER(nbpasmx=1000000)
31	REAL temps(nbpasmx)
32	INTEGER unitlec
33	INTEGER i,j,l,jj
34	REAL constR
35
36	c Declarations NCDF:
37	c -----------------
38	CHARACTER*100 varname
39	INTEGER ierr,nid,nvarid,dimid
40	LOGICAL nc
41	INTEGER start_ps(3),start_temp(4),start_co2ice(3)
42	INTEGER count_ps(3),count_temp(4),count_co2ice(3)
43
44	c declarations pour les points viking:
45	c ------------------------------------
46	INTEGER ivik(2),jvik(2),ifile(2),iv
47
48	REAL, PARAMETER :: lonvik1 = -47.95
49	REAL, PARAMETER :: latvik1 = 22.27
50	REAL, PARAMETER :: lonvik2 = 134.29
51	REAL, PARAMETER :: latvik2 = 47.67
52
53	REAL, PARAMETER :: phivik1 = -3637
54	REAL, PARAMETER :: phivik2 = -4505
55
56
57	REAL lonvik(2),latvik(2),phivik(2),phisim(2)
58	REAL unanj
59
60	c variables meteo:
61	c ----------------
62	REAL vnat(iip1,jjm,llm),unat(iip1,jjp1,llm)
63	REAL t(iip1,jjp1,llm),ps(iip1,jjp1),pstot, phis(iip1,jjp1)
64	REAL co2ice(iip1,jjp1), captotN,captotS
65	real t7(iip1,jjp1) ! temperature in 7th atmospheric layer
66
67	REAL zp1,zp2,zp2_sm,zu,zv,zw(0:1,0:1,2),zalpha,zbeta
68
69	LOGICAL firstcal
70	INTEGER*4 day0
71
72	REAL ziceco2(iip1,jjp1)
73	REAL day,zt,sollong,sol,dayw,dayw_ls
74	REAL airtot1,gh
75
76	INTEGER ii,iyear,kyear
77
78	CHARACTER*2 chr2
79
80
81	c declarations de l'interface avec mywrite:
82	c -----------------------------------------
83
84	CHARACTER file*80
85	CHARACTER pathchmp80,pathsor80,nomfich*80
86
87	INTEGER Time_unit
88
89
90	c externe:
91	c --------
92
93	EXTERNAL iniconst,inigeom,covcont,mywrite
94	EXTERNAL exner,pbar
95	EXTERNAL coordij,moy2
96	EXTERNAL SSUM
97	REAL SSUM
98
99
100	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
101
102	c-----------------------------------------------------------------------
103	c initialisations:
104	c-----------------------------------------------------------------------
105
106	chr2="0"
107	unanj=669.
108	print*,'WARNING!!!',unanj,'Jours/an'
109	nc=.true.
110
111	phivik(1) = phivik1
112	phivik(2) = phivik2
113
114	print *, 'COORDVIKIIIN', latvik, lonvik
115	print*, 'LES PHIVIK', phivik
116
117
118
119
120
121	WRITE(,) 'Chemin des fichiers histoires'
122	READ (*,'(a)') pathchmp
123	WRITE(,) 'Chemin des fichiers sorties'
124	READ (*,'(a)') pathsor
125
126	WRITE(,) 'Fichiers de sortie en sol (1)
127	&,en ls (2) ,les deux (3)'
128	READ (,) Time_unit
129
130
131	write (,)'>>>>>>>>>>>>>>>>', phivik,g
132	DO iv=1,2
133	phivik(iv)=phivik(iv)*3.73
134	END DO
135
136	c-----------------------------------------------------------------------
137	c ouverture des fichiers xgraph:
138	c-----------------------------------------------------------------------
139	ifile(1)=12
140	ifile(2)=13
141	kyear=-1
142	unitlec=11
143
144
145	print*,'Entrer un fichier NC (sans le .nc)'
146	READ(5,'(a)',err=9999) nomfich
147
148
149	c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
150	c grande boucle sur les fichiers histoire:
151	c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
152
153	firstcal=.true.
154	DO WHILE(len_trim(nomfich).GT.0.AND.len_trim(nomfich).LT.50)
155	PRINT *,'>>> nomfich : ',trim(nomfich)
156
157	c----------------------------------------------------------------------
158	c Ouverture des fichiers histoire:
159	c----------------------------------------------------------------------
160
161	file=pathchmp(1:len_trim(pathchmp))//'/'//
162	s nomfich(1:len_trim(nomfich))
163	PRINT*,'file.nc: ', file(1:len_trim(file))//'.nc'
164	PRINT*,'timestep ',dtvr
165
166	IF(nc) THEN
167	ierr= NF_OPEN(file(1:len_trim(file))//'.nc',NF_NOWRITE,nid)
168	ELSE
169	PRINT*,'Ouverture binaire ',file
170	OPEN(unitlec,file=file,status='old',form='unformatted',
171	. iostat=ierr)
172	ENDIF
173
174	c----------------------------------------------------------------------
175	c initialisation de la physique:
176	c----------------------------------------------------------------------
177
178	CALL readhead_NC(file(1:len_trim(file))//'.nc',day0,phis,constR)
179
180	WRITE (,) 'day0 = ' , day0
181
182	CALL conf_gcm( 99, .TRUE. )
183	CALL iniconst
184	CALL inigeom
185
186	c----------------------------------------------------------------------
187	c Lecture temps :
188	c----------------------------------------------------------------------
189
190
191	ierr= NF_INQ_DIMID (nid,"Time",dimid)
192	IF (ierr.NE.NF_NOERR) THEN
193	PRINT*, 'xvik: Le champ <Time> est absent'
194	CALL abort
195	ENDIF
196
197	ierr= NF_INQ_DIMLEN (nid,dimid,nbpas)
198
199	ierr = NF_INQ_VARID (nid, "Time", nvarid)
200	#ifdef NC_DOUBLE
201	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, temps)
202	#else
203	ierr = NF_GET_VAR_REAL(nid, nvarid, temps)
204	#endif
205	IF (ierr.NE.NF_NOERR) THEN
206	PRINT*, 'xvik: Lecture echouee pour <Time>'
207	CALL abort
208	ENDIF
209
210	PRINT*,'temps(1:10)',(temps(itau),itau=1,10)
211
212
213	c-----------------------------------------------------------------------
214	c coordonnees des point Viking:
215	c --------------------------------------------------------------------
216
217	lonvik(1) = lonvik1 * pi/180.
218	latvik(1) = latvik1 * pi/180.
219	lonvik(2) = lonvik2 * pi/180.
220	latvik(2) = latvik2 * pi/180.
221
222
223	c----------------------------------------------------------------------
224	c ponderations pour les 4 points autour de Viking
225	c----------------------------------------------------------------------
226
227
228	DO iv=1,2
229	! locate index of GCM grid points near VL
230	do i=1,iim
231	! we know longitudes are ordered -180...180
232	if ((lonvik(iv).ge.rlonu(i)).and.
233	& (lonvik(iv).le.rlonu(i+1))) then
234	ivik(iv)=i
235	exit
236	endif
237	enddo
238	do j=1,jjm-1
239	!we know tha latitudes are ordered 90...-90
240	if ((latvik(iv).le.rlatv(j)).and.
241	& (latvik(iv).ge.rlatv(j+1))) then
242	jvik(iv)=j
243	exit
244	endif
245	enddo
246	zalpha=(lonvik(iv)-rlonu(ivik(iv)))/
247	s (rlonu(ivik(iv)+1)-rlonu(ivik(iv)))
248	zbeta=(latvik(iv)-rlatv(jvik(iv)))/
249	s (rlatv(jvik(iv)+1)-rlatv(jvik(iv)))
250	zw(0,0,iv)=(1.-zalpha)*(1.-zbeta)
251	zw(1,0,iv)=zalpha*(1.-zbeta)
252	zw(0,1,iv)=(1.-zalpha)*zbeta
253	zw(1,1,iv)=zalpha*zbeta
254	ENDDO
255
256	c----------------------------------------------------------------------
257	c altitude reelle et modele aux points Viking
258	c----------------------------------------------------------------------
259
260
261	DO iv=1,2
262	phisim(iv)=0.
263	DO jj=0,1
264	j=jvik(iv)+jj
265	DO ii=0,1
266	i=ivik(iv)+ii
267	phisim(iv)=phisim(iv)+zw(ii,jj,iv)*phis(i,j)
268	ENDDO
269	ENDDO
270	ENDDO
271	PRINT*,'relief aux points Viking pour les sorties:',phivik
272
273
274	c----------------------------------------------------------------------
275	c lectures des etats:
276	c -------------------------------------------------------------------
277
278	airtot1=1./(SSUM(ip1jmp1,aire,1)-SSUM(jjp1,aire,iip1))
279
280	c======================================================================
281	c debut de la boucle sur les etats dans un fichier histoire:
282	c======================================================================
283
284
285	count_ps=(/iip1,jjp1,1/)
286	count_co2ice=(/iip1,jjp1,1/)
287	count_temp=(/iip1,jjp1,llm,1/)
288
289	DO itau=1,nbpas
290
291	start_ps=(/1,1,itau/)
292	start_co2ice=(/1,1,itau/)
293	start_temp=(/1,1,1,itau/)
294
295	c----------------------------------------------------------------------
296	c lecture drs des champs:
297	c----------------------------------------------------------------------
298
299
300	ccccccccc LECTURE Ps ccccccccccccccccccccccccccc
301
302
303	ierr = NF_INQ_VARID (nid, "ps", nvarid)
304	#ifdef NC_DOUBLE
305	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start_ps,count_ps, ps)
306	#else
307	ierr = NF_GET_VARA_REAL(nid, nvarid,start_ps,count_ps, ps)
308	#endif
309	IF (ierr.NE.NF_NOERR) THEN
310	PRINT*, 'xvik: Lecture echouee pour <ps>'
311	CALL abort
312	ENDIF
313
314	PRINT*,'ps',ps(iip1/2,jjp1/2)
315
316	ccccccccc LECTURE Temperature ccccccccccccccccccccccccccc
317
318
319	ierr = NF_INQ_VARID (nid, "temp", nvarid)
320	#ifdef NC_DOUBLE
321	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start_temp,count_temp, t)
322	#else
323	ierr = NF_GET_VARA_REAL(nid,nvarid,start_temp,count_temp, t)
324	#endif
325	IF (ierr.NE.NF_NOERR) THEN
326	PRINT*, 'xvik: Lecture echouee pour <temp>'
327	! Ehouarn: proceed anyways
328	! CALL abort
329	write(,)'--> Setting temperature to zero !!!'
330	t(1:iip1,1:jjp1,1:llm)=0.0
331	write(,)'--> looking for temp7 (temp in 7th layer)'
332	ierr=NF_INQ_VARID(nid,"temp7", nvarid)
333	if (ierr.eq.NF_NOERR) then
334	write(,) " OK, found temp7 variable"
335	#ifdef NC_DOUBLE
336	ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start_ps,count_ps,t7)
337	#else
338	ierr=NF_GET_VARA_REAL(nid,nvarid,start_ps,count_ps,t7)
339	#endif
340	if (ierr.ne.NF_NOERR) then
341	write(,)'xvik: failed loading temp7 !'
342	stop
343	endif
344	else ! no 'temp7' variable
345	write(,)' No temp7 variable either !'
346	write(,)' Will have to to without ...'
347	t7(1:iip1,1:jjp1)=0.0
348	endif
349	ELSE ! t() was successfully loaded, copy 7th layer to t7()
350	t7(1:iip1,1:jjp1)=t(1:iip1,1:jjp1,7)
351	ENDIF
352
353
354
355	ccccccccc LECTURE co2ice ccccccccccccccccccccccccccc
356
357
358	ierr = NF_INQ_VARID (nid, "co2ice", nvarid)
359	#ifdef NC_DOUBLE
360	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start_co2ice,
361	& count_co2ice, co2ice)
362	#else
363	ierr = NF_GET_VARA_REAL(nid, nvarid,start_co2ice,
364	& count_co2ice, co2ice)
365	#endif
366	IF (ierr.NE.NF_NOERR) THEN
367	PRINT*, 'xvik: Lecture echouee pour <co2ice>'
368	CALL abort
369	ENDIF
370
371	c----------------------------------------------------------------------
372	c Gestion du temps
373	c ---------------------------------------------------------------------
374
375	day=temps(itau)
376	PRINT*,'day ',day
377	sol=day+day0
378	iyear=sol/unanj
379	WRITE (,) 'iyear',iyear
380	sol=sol-iyear*unanj
381
382	c----------------------------------------------------------------------
383	c Ouverture / fermeture des fichiers
384	c ---------------------------------------------------------------------
385
386	IF (iyear.NE.kyear) THEN
387	WRITE(chr2(1:1),'(i1)') iyear+1
388	WRITE (,) 'iyear bis',iyear
389	WRITE (,) 'chr2'
390	WRITE (,) chr2
391	IF(iyear.GE.9) WRITE(chr2,'(i2)') iyear+1
392	kyear=iyear
393	DO ii=1,2
394	CLOSE(10+ifile(ii))
395	CLOSE(2+ifile(ii))
396	CLOSE(4+ifile(ii))
397	CLOSE(6+ifile(ii))
398	CLOSE(8+ifile(ii))
399	CLOSE(16+ifile(ii))
400	CLOSE(12+ifile(ii))
401	CLOSE(14+ifile(ii))
402	CLOSE(97)
403	CLOSE(98)
404	ENDDO
405	CLOSE(5+ifile(1))
406	OPEN(ifile(1)+10,file='xpsol1'//chr2,form='formatted')
407	OPEN(ifile(2)+10,file='xpsol2'//chr2,form='formatted')
408	OPEN(97,file='xprestot'//chr2,form='formatted')
409
410	ENDIF
411
412	dayw = sol
413	call sol2ls(sol,sollong)
414	dayw_ls = sollong
415
416
417
418	c----------------------------------------------------------------------
419	c Calcul de la moyenne de pression planetaire
420	c ---------------------------------------------------------------------
421
422
423	pstot=0.
424	captotS=0.
425	captotN=0.
426	DO j=1,jjp1
427	DO i=1,iim
428	pstot=pstot+aire(i,j)*ps(i,j)
429	ENDDO
430	ENDDO
431
432	DO j=1,jjp1/2
433	DO i=1,iim
434	captotN = captotN +aire(i,j)*co2ice(i,j)
435	ENDDO
436	ENDDO
437	DO j=jjp1/2+1, jjp1
438	DO i=1,iim
439	captotS = captotS +aire(i,j)*co2ice(i,j)
440	ENDDO
441	ENDDO
442
443
444	c --------------Ecriture fichier sortie xprestot-----------------------
445	c Sol ou ls ou les deux
446	c Ps_moy_planetaire (Pa)
447	c Pequivalente de glace de CO2 au Nord (si entierement sublimee) (Pa)
448	c Pequivalente de glace de CO2 au Sud (si entierement sublimee) (Pa)
449
450
451	IF(Time_unit == 1) THEN
452	WRITE(97,'(4e16.6)') dayw,pstot*airtot1
453	& , captotNgairtot1, captotSgairtot1
454
455	ELSEIF (Time_unit == 2) THEN
456	WRITE(97,'(4e16.6)') dayw_ls,pstot*airtot1
457	& , captotNgairtot1, captotSgairtot1
458
459	ELSE
460	WRITE(97,'(5e16.6)') dayw,dayw_ls,pstot*airtot1
461	& , captotNgairtot1,captotSgairtot1
462
463
464	ENDIF
465
466	c----------------------------------------------------------------------
467	c boucle sur les sites vikings:
468	c----------------------------------------------------------------------
469
470	c----------------------------------------------------------------------
471	c interpolation de la temperature dans la 7eme couche, de la pression
472	c de surface et des vents aux points viking.
473	c----------------------------------------------------------------------
474
475	IF(.NOT.firstcal) THEN
476
477	DO iv=1,2
478
479	zp1=0.
480	zp2=0.
481	zp2_sm=0.
482	zt=0.
483
484	DO jj=0,1
485
486	j=jvik(iv)+jj
487
488	DO ii=0,1
489
490	i=ivik(iv)+ii
491	zt=zt+zw(ii,jj,iv)*t7(i,j)
492	zp1=zp1+zw(ii,jj,iv)*log(ps(i,j)) ! interpolate in log(P)
493	WRITE (,) 'ps autour iv',ps(i,j),iv
494
495	ENDDO
496	ENDDO
497
498	zp1=exp(zp1) ! because of the bilinear interpolation in log(P)
499	WRITE (,) 'constR ',constR
500	WRITE (,) 'zt ',zt
501	gh=constR*zt
502
503	c----------------------------------------------------------------------
504	c pression au sol extrapolee a partir de la temp. 7eme couche
505	c----------------------------------------------------------------------
506
507	if (gh.eq.0) then ! if we don't have temperature values
508	! assume a scale height of 10km
509	zp2=zp1exp(-(phivik(iv)-phisim(iv))/(3.731.e4))
510	else
511	zp2=zp1*exp(-(phivik(iv)-phisim(iv))/gh)
512	endif
513
514	WRITE (,) 'iv,pstot,zp2, zp1, phivik(iv),phisim(iv),gh'
515	WRITE (,) iv,pstot*airtot1,zp2,zp1,phivik(iv),phisim(iv),gh
516
517
518	c ------Ecriture 2 fichiers (1 pour Vl1, 1 pour VL2) sortie xpsol ------
519	c Sol ou ls ou les deux
520	c Ps site VLi (i=1,2) a l'altitude GCM (Pa)
521	c Ps site VLi (i=1,2) a l'altitude exacte (interpolee) (Pa)
522
523	IF(Time_unit == 1) THEN
524	WRITE(ifile(iv)+10,'(3e15.5)') dayw,zp2,zp1
525	ELSEIF (Time_unit == 2) THEN
526	WRITE(ifile(iv)+10,'(3e15.5)') dayw_ls,zp2,zp1
527	ELSE
528	WRITE(ifile(iv)+10,'(4e15.5)') dayw,dayw_ls,zp2,zp1
529	ENDIF
530
531	ENDDO
532
533	ENDIF
534
535	firstcal=.false.
536
537
538	c======================================================================
539	c Fin de la boucle sur les etats du fichier histoire:
540	c======================================================================
541
542	ENDDO
543
544	ierr= NF_CLOSE(nid)
545
546	PRINT*,'Fin du fichier',nomfich
547	print*,'Entrer un nouveau fichier NC
548	&(sans le .nc) ou return pour finir'
549	READ(5,'(a)',err=9999) nomfich
550
551
552	c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
553	c Fin de la boucle sur les fichiers histoire:
554	c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
555
556	ENDDO
557
558	PRINT,'relief du point V1',.001phis(ivik(1),jvik(1))/g
559	PRINT,'relief du point V2',.001phis(ivik(2),jvik(2))/g
560	DO iv=1,2
561	PRINT*,'Viking',iv,' i=',ivik(iv),'j =',jvik(iv)
562	WRITE(6,7777)
563	s (rlonv(i)*180./pi,i=ivik(iv)-1,ivik(iv)+2)
564	print*
565	DO j=jvik(iv)-1,jvik(iv)+2
566	WRITE(6,'(f8.1,10x,5f7.1)')
567	s rlatu(j)180./pi,(phis(i,j)/(g1000.),i=ivik(iv)-1,ivik(iv)+2)
568	ENDDO
569	print*
570	print*,'zw'
571	write(6,'(2(2f10.4/))') ((zw(ii,jj,iv),ii=0,1),jj=0,1)
572	print*,'altitude interpolee (km) ',phisim(iv)/1000./g
573	ENDDO
574	PRINT*,'R=',r
575	9999 PRINT*,'Fin '
576
577	7777 FORMAT ('latitude/longitude',4f7.1)
578
579
580
581	END
582
583	subroutine sol2ls(sol,Ls)
584	!==============================================================================
585	! Purpose:
586	! Convert a date/time, given in sol (martian day),
587	! into solar longitude date/time, in Ls (in degrees),
588	! where sol=0 is (by definition) the northern hemisphere
589	! spring equinox (where Ls=0).
590	!==============================================================================
591	! Notes:
592	! Even though "Ls" is cyclic, if "sol" is greater than N (martian) year,
593	! "Ls" will be increased by N*360
594	! Won't work as expected if sol is negative (then again,
595	! why would that ever happen?)
596	!==============================================================================
597
598	implicit none
599
600	!==============================================================================
601	! Arguments:
602	!==============================================================================
603	real,intent(in) :: sol
604	real,intent(out) :: Ls
605
606	!==============================================================================
607	! Local variables:
608	!==============================================================================
609	real year_day,peri_day,timeperi,e_elips,twopi,degrad
610	data year_day /669./ ! # of sols in a martian year
611	data peri_day /485.0/
612	data timeperi /1.9082314/
613	data e_elips /0.093358/
614	data twopi /6.2831853/ ! 2.*pi
615	data degrad /57.2957795/ ! pi/180
616
617	real zanom,xref,zx0,zdx,zteta,zz
618
619	integer count_years
620	integer iter
621
622	!==============================================================================
623	! 1. Compute Ls
624	!==============================================================================
625
626	zz=(sol-peri_day)/year_day
627	zanom=twopi*(zz-nint(zz))
628	xref=abs(zanom)
629
630	! The equation zx0 - e * sin (zx0) = xref, solved by Newton
631	zx0=xref+e_elips*sin(xref)
632	do iter=1,20 ! typically, 2 or 3 iterations are enough
633	zdx=-(zx0-e_elipssin(zx0)-xref)/(1.-e_elipscos(zx0))
634	zx0=zx0+zdx
635	if(abs(zdx).le.(1.e-7)) then
636	! write(,)'iter:',iter,' \|zdx\|:',abs(zdx)
637	exit
638	endif
639	enddo
640
641	if(zanom.lt.0.) zx0=-zx0
642
643	zteta=2.atan(sqrt((1.+e_elips)/(1.-e_elips))tan(zx0/2.))
644	Ls=zteta-timeperi
645
646	if(Ls.lt.0.) then
647	Ls=Ls+twopi
648	else
649	if(Ls.gt.twopi) then
650	Ls=Ls-twopi
651	endif
652	endif
653
654	Ls=degrad*Ls
655	! Ls is now in degrees
656
657	!==============================================================================
658	! 1. Account for (eventual) years included in input date/time sol
659	!==============================================================================
660
661	count_years=0 ! initialize
662	zz=sol ! use "zz" to store (and work on) the value of sol
663	do while (zz.ge.year_day)
664	count_years=count_years+1
665	zz=zz-year_day
666	enddo
667
668	! Add 360 degrees to Ls for every year
669	if (count_years.ne.0) then
670	Ls=Ls+360.*count_years
671	endif
672
673
674	end subroutine sol2ls

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

Download in other formats:

Original Format