Context Navigation

xvik.F @ 3094

Last change on this file since 3094 was 2844, checked in by tpierron, 2 years ago

Mars GCM :
Update xvik program. The program now interpolates surface pressure at 3 locations : Viking Lander 1, Viking Lander 2 and Insight.
Also add as outputs :

Diurnal average
Harmonics fit (Fourier series) of the diurnal average (6 harmonics)

Columns of each output are described at the top of each file. (Fourier coefficients are also given for harmonics fit)
Total CO2 inventory is also caluclated now by xvik.
TP

File size: 38.5 KB

Rev	Line
[38]	1	PROGRAM xvik
[1403]	2
	3	USE filtreg_mod, ONLY: inifilr
[1422]	4	USE comconst_mod, ONLY: dtvr,g,r,pi
[2824]	5	USE comvert_mod, ONLY: pa,preff
[2325]	6
[38]	7	IMPLICIT NONE
[2325]	8
	9
[38]	10	c=======================================================================
	11	c
[2844]	12	c Pressure at Insight and Viking sites
[38]	13	c
	14	c=======================================================================
[2325]	15
	16
[38]	17	c-----------------------------------------------------------------------
	18	c declarations:
[2325]	19	c-----------------------------------------------------------------------
[38]	20
	21
[1977]	22	include "dimensions.h"
	23	include "paramet.h"
	24	include "comdissip.h"
	25	include "comgeom2.h"
	26	include "netcdf.inc"
[38]	27
	28
[2325]	29	INTEGER itau,nbpas,nbpasmx
[38]	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	INTEGER start_ps(3),start_temp(4),start_co2ice(3)
	41	INTEGER count_ps(3),count_temp(4),count_co2ice(3)
	42
	43	c declarations pour les points viking:
	44	c ------------------------------------
[2844]	45	INTEGER isite(3),jsite(3),ifile(3),iv
[2325]	46
	47	REAL, PARAMETER :: lonvik1 = -47.95
	48	REAL, PARAMETER :: latvik1 = 22.27
	49	REAL, PARAMETER :: lonvik2 = 134.29
	50	REAL, PARAMETER :: latvik2 = 47.67
[2844]	51	REAL, PARAMETER :: loninst = 135.62
	52	REAL, PARAMETER :: latinst = 4.502
	53
[2325]	54	REAL, PARAMETER :: phivik1 = -3637
	55	REAL, PARAMETER :: phivik2 = -4505
[2844]	56	REAL, PARAMETER :: phiinst = -2614
[2325]	57
	58
[2844]	59	REAL lonsite(3),latsite(3),phisite(3),phisim(3)
[38]	60	REAL unanj
	61
	62	c variables meteo:
	63	c ----------------
	64	REAL vnat(iip1,jjm,llm),unat(iip1,jjp1,llm)
	65	REAL t(iip1,jjp1,llm),ps(iip1,jjp1),pstot, phis(iip1,jjp1)
	66	REAL co2ice(iip1,jjp1), captotN,captotS
	67	real t7(iip1,jjp1) ! temperature in 7th atmospheric layer
	68
[2844]	69	REAL zp1,zp2,zp2_sm,zu,zv,zw(0:1,0:1,3),zalpha,zbeta
[38]	70
[2325]	71	LOGICAL firstcal
[38]	72	INTEGER*4 day0
	73
	74	REAL ziceco2(iip1,jjp1)
[2325]	75	REAL day,zt,sollong,sol,dayw,dayw_ls
[38]	76	REAL airtot1,gh
	77
	78	INTEGER ii,iyear,kyear
	79
	80	CHARACTER*2 chr2
	81
	82
	83	c declarations de l'interface avec mywrite:
	84	c -----------------------------------------
	85
	86	CHARACTER file*80
	87	CHARACTER pathchmp80,pathsor80,nomfich*80
[2325]	88
	89	INTEGER Time_unit
	90
[2844]	91	REAL ls2sol
	92
[38]	93
	94	c externe:
	95	c --------
	96
	97	EXTERNAL iniconst,inigeom,covcont,mywrite
[1403]	98	EXTERNAL exner,pbar
[2325]	99	EXTERNAL coordij,moy2
[38]	100	EXTERNAL SSUM
	101	REAL SSUM
[2844]	102
	103	c diurnam:
	104	c --------
	105	integer di,di_prev
	106	integer k
	107	real ps_gcm_diurnal, ps_diurnal
	108	integer compt_diurn
	109
	110	c harmonics:
	111	c --------
	112
	113	integer, parameter :: nbmax = 999999
	114	integer ndata
	115	real ls_harm
	116	real ls_tab(nbmax)
	117	real ps_diurnal_tab(nbmax),ps_gcm_diurnal_tab(nbmax)
	118	real a_tab(nbmax),b_tab(nbmax)
	119	real a_tab_gcm(nbmax),b_tab_gcm(nbmax)
	120	real a,b
	121	real ps_harm, ps_gcm_harm
	122	integer, parameter :: n_harmo = 6
[2325]	123
[2844]	124
[38]	125	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	126
	127	c-----------------------------------------------------------------------
	128	c initialisations:
[2325]	129	c-----------------------------------------------------------------------
[2824]	130	pi=4.*atan(1.)
	131	pa=20
	132	preff=610.
[38]	133
[1977]	134	chr2="0"
[2824]	135	iyear=0
[1977]	136	unanj=669.
[2824]	137	print*,'WARNING!!! Assuming',unanj,'sols/year'
[2325]	138
[2824]	139	c-----------------------------------------------------------------------
[2844]	140	c Viking Lander and Insight coordinates:
[2824]	141	c --------------------------------------------------------------------
	142
[2844]	143	lonsite(1) = lonvik1
	144	latsite(1) = latvik1
	145	lonsite(2) = lonvik2
	146	latsite(2) = latvik2
	147	lonsite(3) = loninst
	148	latsite(3) = latinst
[2824]	149
[2844]	150	phisite(1) = phivik1
	151	phisite(2) = phivik2
	152	phisite(3) = phiinst
	153
	154	WRITE(,) 'Viking1 coordinates:'
	155	WRITE(,) 'latvik1:',latvik1,' lonvik1:',lonvik1
	156	WRITE(,) 'Phivik1:', phivik1
[2325]	157
[2844]	158	WRITE(,) 'Viking2 coordinates:'
	159	WRITE(,) 'latvik2:',latvik2,' lonvik2:',lonvik2
	160	WRITE(,) 'Phivik2:', phivik2
[2325]	161
[2844]	162	WRITE(,) 'Insight coordinates:'
	163	WRITE(,) 'latinst:',latinst,' loninst:',loninst
	164	WRITE(,) 'Phiinst:', phiinst
	165
[2824]	166	! convert coordinates to radians
[2844]	167	lonsite(1) = lonvik1 * pi/180.
	168	latsite(1) = latvik1 * pi/180.
	169	lonsite(2) = lonvik2 * pi/180.
	170	latsite(2) = latvik2 * pi/180.
	171	lonsite(3) = loninst * pi/180.
	172	latsite(3) = latinst * pi/180.
[2325]	173
[2844]	174
	175
[2824]	176	WRITE(,) 'Path to the diagfi files directory'
[2325]	177	READ (*,'(a)') pathchmp
[2824]	178	WRITE(,) 'Path to the dir for outputs'
[2325]	179	READ (*,'(a)') pathsor
	180
[2824]	181	WRITE(,) 'Output file time axis in sol (1) '//
	182	&'in ls (2) ,or both (3)'
[2325]	183	READ (,) Time_unit
	184
	185
[2844]	186	write (,)'>>>>>>>>>>>>>>>>', phisite,g
	187	DO iv=1,3
	188	phisite(iv)=phisite(iv)*3.73
[38]	189	END DO
	190
	191	c-----------------------------------------------------------------------
[2824]	192	c output files:
[38]	193	c-----------------------------------------------------------------------
	194	ifile(1)=12
	195	ifile(2)=13
[2844]	196	ifile(3)=14
	197
[38]	198	kyear=-1
	199	unitlec=11
[2325]	200
	201
[2824]	202	print*,'diagfi file name (without trailing .nc)'
[2325]	203	READ(5,'(a)',err=9999) nomfich
	204
[38]	205
	206	c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[2824]	207	c loop on the diagfi files:
[38]	208	c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	209
	210	firstcal=.true.
[1403]	211	DO WHILE(len_trim(nomfich).GT.0.AND.len_trim(nomfich).LT.50)
[38]	212	PRINT *,'>>> nomfich : ',trim(nomfich)
	213
[2325]	214	c----------------------------------------------------------------------
[2844]	215	c Open diagfi files :
[2325]	216	c----------------------------------------------------------------------
[38]	217
[1403]	218	file=pathchmp(1:len_trim(pathchmp))//'/'//
	219	s nomfich(1:len_trim(nomfich))
	220	PRINT*,'file.nc: ', file(1:len_trim(file))//'.nc'
[38]	221	PRINT*,'timestep ',dtvr
	222
[1403]	223	ierr= NF_OPEN(file(1:len_trim(file))//'.nc',NF_NOWRITE,nid)
[38]	224
	225	c----------------------------------------------------------------------
[2824]	226	c initialise physics:
[2325]	227	c----------------------------------------------------------------------
[38]	228
[1403]	229	CALL readhead_NC(file(1:len_trim(file))//'.nc',day0,phis,constR)
[38]	230
	231	WRITE (,) 'day0 = ' , day0
	232
[1977]	233	CALL conf_gcm( 99, .TRUE. )
[38]	234	CALL iniconst
	235	CALL inigeom
	236
[2325]	237	c----------------------------------------------------------------------
[2844]	238	c Read time :
[2325]	239	c----------------------------------------------------------------------
[38]	240
[2325]	241
[38]	242	ierr= NF_INQ_DIMID (nid,"Time",dimid)
	243	IF (ierr.NE.NF_NOERR) THEN
	244	PRINT*, 'xvik: Le champ <Time> est absent'
	245	CALL abort
	246	ENDIF
	247
	248	ierr= NF_INQ_DIMLEN (nid,dimid,nbpas)
	249
	250	ierr = NF_INQ_VARID (nid, "Time", nvarid)
	251	#ifdef NC_DOUBLE
	252	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, temps)
	253	#else
	254	ierr = NF_GET_VAR_REAL(nid, nvarid, temps)
	255	#endif
	256	IF (ierr.NE.NF_NOERR) THEN
	257	PRINT*, 'xvik: Lecture echouee pour <Time>'
	258	CALL abort
	259	ENDIF
	260
[1977]	261	PRINT*,'temps(1:10)',(temps(itau),itau=1,10)
[2325]	262
	263
	264
[2844]	265	c------------------------------------------------------
	266	c Weights for 4 near points at Viking and Insight
	267	c------------------------------------------------------
[38]	268
[2844]	269	DO iv=1,3
[1977]	270	! locate index of GCM grid points near VL
[2844]	271	do i=1,iim
[1977]	272	! we know longitudes are ordered -180...180
[2844]	273	write(,) i, lonsite(iv),rlonu(i),rlonu(i+1)
	274	if ((lonsite(iv).ge.rlonu(i)).and.
	275	& (lonsite(iv).le.rlonu(i+1))) then
	276	isite(iv)=i
[1977]	277	exit
	278	endif
	279	enddo
[2844]	280
[1977]	281	do j=1,jjm-1
	282	!we know tha latitudes are ordered 90...-90
[2844]	283	if ((latsite(iv).le.rlatv(j)).and.
	284	& (latsite(iv).ge.rlatv(j+1))) then
	285	jsite(iv)=j
[1977]	286	exit
	287	endif
	288	enddo
[2844]	289	zalpha=(lonsite(iv)-rlonu(isite(iv)))/
	290	s (rlonu(isite(iv)+1)-rlonu(isite(iv)))
	291	zbeta=(latsite(iv)-rlatv(jsite(iv)))/
	292	s (rlatv(jsite(iv)+1)-rlatv(jsite(iv)))
[38]	293	zw(0,0,iv)=(1.-zalpha)*(1.-zbeta)
	294	zw(1,0,iv)=zalpha*(1.-zbeta)
	295	zw(0,1,iv)=(1.-zalpha)*zbeta
	296	zw(1,1,iv)=zalpha*zbeta
	297	ENDDO
	298
[2325]	299	c----------------------------------------------------------------------
[2844]	300	c true and model altitude at Viking and Insight locations
[2325]	301	c----------------------------------------------------------------------
	302
	303
[2844]	304	DO iv=1,3
[38]	305	phisim(iv)=0.
	306	DO jj=0,1
[2844]	307	j=jsite(iv)+jj
[38]	308	DO ii=0,1
[2844]	309	i=isite(iv)+ii
[38]	310	phisim(iv)=phisim(iv)+zw(ii,jj,iv)*phis(i,j)
	311	ENDDO
	312	ENDDO
	313	ENDDO
[2844]	314	PRINT*,'phisite at Viking locations for outputs:',phisite
[2325]	315
[38]	316
	317	c----------------------------------------------------------------------
[2824]	318	c read variables:
[2325]	319	c -------------------------------------------------------------------
[38]	320
	321	airtot1=1./(SSUM(ip1jmp1,aire,1)-SSUM(jjp1,aire,iip1))
	322
	323	c======================================================================
[2844]	324	c Begin the loop on states in inputs files :
[38]	325	c======================================================================
[2325]	326
	327
[38]	328	count_ps=(/iip1,jjp1,1/)
	329	count_co2ice=(/iip1,jjp1,1/)
	330	count_temp=(/iip1,jjp1,llm,1/)
	331
	332	DO itau=1,nbpas
	333
	334	start_ps=(/1,1,itau/)
	335	start_co2ice=(/1,1,itau/)
	336	start_temp=(/1,1,1,itau/)
[2325]	337
	338	c----------------------------------------------------------------------
[2824]	339	c read fields:
[2325]	340	c----------------------------------------------------------------------
[38]	341
[2325]	342
[2824]	343	ccccccccc Load Ps ccccccccccccccccccccccccccc
[2325]	344
	345
[38]	346	ierr = NF_INQ_VARID (nid, "ps", nvarid)
	347	#ifdef NC_DOUBLE
	348	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start_ps,count_ps, ps)
	349	#else
	350	ierr = NF_GET_VARA_REAL(nid, nvarid,start_ps,count_ps, ps)
	351	#endif
	352	IF (ierr.NE.NF_NOERR) THEN
	353	PRINT*, 'xvik: Lecture echouee pour <ps>'
	354	CALL abort
	355	ENDIF
	356
	357	PRINT*,'ps',ps(iip1/2,jjp1/2)
	358
[2824]	359	ccccccccc Load Temperature ccccccccccccccccccccccccccc
[2325]	360
	361
[38]	362	ierr = NF_INQ_VARID (nid, "temp", nvarid)
	363	#ifdef NC_DOUBLE
	364	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start_temp,count_temp, t)
	365	#else
	366	ierr = NF_GET_VARA_REAL(nid,nvarid,start_temp,count_temp, t)
	367	#endif
	368	IF (ierr.NE.NF_NOERR) THEN
	369	PRINT*, 'xvik: Lecture echouee pour <temp>'
	370	! Ehouarn: proceed anyways
	371	! CALL abort
	372	write(,)'--> Setting temperature to zero !!!'
	373	t(1:iip1,1:jjp1,1:llm)=0.0
	374	write(,)'--> looking for temp7 (temp in 7th layer)'
	375	ierr=NF_INQ_VARID(nid,"temp7", nvarid)
	376	if (ierr.eq.NF_NOERR) then
	377	write(,) " OK, found temp7 variable"
	378	#ifdef NC_DOUBLE
	379	ierr=NF_GET_VARA_DOUBLE(nid,nvarid,start_ps,count_ps,t7)
	380	#else
	381	ierr=NF_GET_VARA_REAL(nid,nvarid,start_ps,count_ps,t7)
	382	#endif
	383	if (ierr.ne.NF_NOERR) then
	384	write(,)'xvik: failed loading temp7 !'
	385	stop
	386	endif
	387	else ! no 'temp7' variable
	388	write(,)' No temp7 variable either !'
	389	write(,)' Will have to to without ...'
	390	t7(1:iip1,1:jjp1)=0.0
	391	endif
	392	ELSE ! t() was successfully loaded, copy 7th layer to t7()
	393	t7(1:iip1,1:jjp1)=t(1:iip1,1:jjp1,7)
	394	ENDIF
	395
	396
[2325]	397
[2824]	398	ccccccccc Load co2ice ccccccccccccccccccccccccccc
[2325]	399
	400
[38]	401	ierr = NF_INQ_VARID (nid, "co2ice", nvarid)
	402	#ifdef NC_DOUBLE
	403	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start_co2ice,
	404	& count_co2ice, co2ice)
	405	#else
	406	ierr = NF_GET_VARA_REAL(nid, nvarid,start_co2ice,
	407	& count_co2ice, co2ice)
	408	#endif
	409	IF (ierr.NE.NF_NOERR) THEN
	410	PRINT*, 'xvik: Lecture echouee pour <co2ice>'
	411	CALL abort
	412	ENDIF
	413
[2325]	414	c----------------------------------------------------------------------
[2824]	415	c Handle calendar
[2325]	416	c ---------------------------------------------------------------------
[38]	417
	418	day=temps(itau)
	419	PRINT*,'day ',day
[2325]	420	sol=day+day0
[2824]	421	do while (sol.gt.unanj)
	422	sol=sol-unanj
	423	enddo
	424	WRITE (,) 'sol: ',sol,' iyear:',iyear
[2325]	425
	426	c----------------------------------------------------------------------
[2824]	427	c Open /close files
[2325]	428	c ---------------------------------------------------------------------
	429
[38]	430	IF (iyear.NE.kyear) THEN
	431	WRITE(chr2(1:1),'(i1)') iyear+1
[2824]	432	WRITE (,) 'iyear bis:',iyear
	433	WRITE (,) 'chr2:',trim(chr2)
[38]	434	IF(iyear.GE.9) WRITE(chr2,'(i2)') iyear+1
	435	kyear=iyear
[2844]	436	DO ii=1,3
[38]	437	CLOSE(10+ifile(ii))
[2844]	438	CLOSE(20+ifile(ii))
[38]	439	CLOSE(2+ifile(ii))
	440	CLOSE(4+ifile(ii))
	441	CLOSE(6+ifile(ii))
	442	CLOSE(8+ifile(ii))
	443	CLOSE(16+ifile(ii))
	444	CLOSE(12+ifile(ii))
	445	CLOSE(14+ifile(ii))
	446	CLOSE(97)
	447	CLOSE(98)
	448	ENDDO
	449	CLOSE(5+ifile(1))
[2844]	450	OPEN(ifile(1)+10,file='ps_VL1_year'//chr2,form='formatted')
	451	OPEN(ifile(2)+10,file='ps_VL2_year'//chr2,form='formatted')
	452	OPEN(ifile(3)+10,file='ps_INS_year'//chr2,form='formatted')
	453	OPEN(97,file='prestot_year'//chr2,form='formatted')
	454
[2325]	455
[2844]	456	c Sol or ls or both
	457	c Planetary mean surface pressure (Pa)
	458	c Equivalent pressure of CO2 ice at North Polar cap (Pa)
	459	c Equivalent pressure of CO2 ice at South Polar cap (Pa)
	460	c Total amount of CO2 on the planet (Pa)
	461
	462	IF (Time_unit == 1) THEN
	463
	464	WRITE(ifile(1)+10,'(a)') '# Sol , Surface Pressure at VL1 at
	465	& true (interpolated) altitude (Pa) ,
	466	& Surface Pressure at VL1 at GCM altitude (Pa)'
	467
	468	WRITE(ifile(2)+10,'(a)') '# Sol , Surface Pressure at VL2 at
	469	& true (interpolated) altitude (Pa) ,
	470	& Surface Pressure at VL2 at GCM altitude (Pa)'
	471
	472	WRITE(ifile(3)+10,'(a)') '# Sol , Surface Pressure at Insight at
	473	& true (interpolated) altitude (Pa) ,
	474	& Surface Pressure at Insight at GCM altitude (Pa)'
	475
	476	WRITE(97,'(a)') '# Sol , Planetary Mean Surface Pressure (Pa) ,
	477	& Equivalent pressure of CO2 ice at North Polar cap (Pa) ,
	478	& Equivalent pressure of CO2 ice at South Polar cap (Pa) ,
	479	& Total amount of CO2 on the planet (Pa)'
	480
	481	ELSEIF (Time_unit == 2) THEN
	482
	483	WRITE(ifile(1)+10,'(a)') '# Ls (deg) , Surface Pressure at VL1 at
	484	& true (interpolated) altitude (Pa) ,
	485	& Surface Pressure at VL1 at GCM altitude (Pa)'
	486
	487	WRITE(ifile(2)+10,'(a)') '# Ls (deg) , Surface Pressure at VL2 at
	488	& true (interpolated) altitude (Pa) ,
	489	& Surface Pressure at VL2 at GCM altitude (Pa)'
	490
	491	WRITE(ifile(3)+10,'(a)') '# Ls (deg) , Surface Pressure at Insight at
	492	& true (interpolated) altitude (Pa) ,
	493	& Surface Pressure at Insight at GCM altitude (Pa)'
	494
	495	WRITE(97,'(a)') '# Ls (deg) , Planetary Mean Surface Pressure (Pa) ,
	496	& Equivalent pressure of CO2 ice at North Polar cap (Pa) ,
	497	& Equivalent pressure of CO2 ice at South Polar cap (Pa) ,
	498	& Total amount of CO2 on the planet (Pa)'
	499
	500	ELSE
	501
	502	WRITE(ifile(1)+10,'(a)') '# Sol , Ls (deg) , Surface Pressure at VL1 at
	503	& true (interpolated) altitude (Pa) ,
	504	& Surface Pressure at VL1 at GCM altitude (Pa)'
	505
	506	WRITE(ifile(2)+10,'(a)') '# Sol , Ls (deg) , Surface Pressure at VL2 at
	507	& true (interpolated) altitude (Pa) ,
	508	& Surface Pressure at VL2 at GCM altitude (Pa)'
	509
	510	WRITE(ifile(3)+10,'(a)') '# Sol , Ls (deg) , Surface Pressure at Insight at
	511	& true (interpolated) altitude (Pa) ,
	512	& Surface Pressure at Insight at GCM altitude (Pa)'
	513
	514	WRITE(97,'(a)') '# Sol , Ls (deg) , Planetary Mean Surface Pressure (Pa) ,
	515	& Equivalent pressure of CO2 ice at North Polar cap (Pa) ,
	516	& Equivalent pressure of CO2 ice at South Polar cap (Pa) ,
	517	& Total amount of CO2 on the planet (Pa)'
	518
	519	ENDIF
	520
	521
[38]	522	ENDIF
	523
[2325]	524	dayw = sol
	525	call sol2ls(sol,sollong)
	526	dayw_ls = sollong
	527
	528
	529
	530	c----------------------------------------------------------------------
[2824]	531	c Compute average planetary pressure
[2325]	532	c ---------------------------------------------------------------------
[38]	533
	534
	535	pstot=0.
	536	captotS=0.
	537	captotN=0.
	538	DO j=1,jjp1
	539	DO i=1,iim
	540	pstot=pstot+aire(i,j)*ps(i,j)
	541	ENDDO
	542	ENDDO
	543
	544	DO j=1,jjp1/2
	545	DO i=1,iim
	546	captotN = captotN +aire(i,j)*co2ice(i,j)
	547	ENDDO
	548	ENDDO
	549	DO j=jjp1/2+1, jjp1
	550	DO i=1,iim
	551	captotS = captotS +aire(i,j)*co2ice(i,j)
	552	ENDDO
	553	ENDDO
[2325]	554
	555
[2844]	556	c --------------Write output file prestot-----------------------
	557	c Sol or ls or both
	558	c Planetary mean surface pressure (Pa)
	559	c Equivalent pressure of CO2 ice at North Polar cap (Pa)
	560	c Equivalent pressure of CO2 ice at South Polar cap (Pa)
	561	c Total amount of CO2 on the planet (Pa)
[2325]	562
	563
	564	IF(Time_unit == 1) THEN
[2844]	565	WRITE(97,'(5e16.6)') dayw,pstot*airtot1
	566	& , captotNgairtot1, captotSgairtot1,
	567	& pstotairtot1+captotNgairtot1+captotSg*airtot1
[2325]	568
	569	ELSEIF (Time_unit == 2) THEN
[2844]	570	WRITE(97,'(5e16.6)') dayw_ls,pstot*airtot1
	571	& , captotNgairtot1, captotSgairtot1,
	572	& pstotairtot1+captotNgairtot1+captotSg*airtot1
[2325]	573
	574	ELSE
[2844]	575	WRITE(97,'(6e16.6)') dayw,dayw_ls,pstot*airtot1
	576	& , captotNgairtot1,captotSgairtot1,
	577	& pstotairtot1+captotNgairtot1+captotSg*airtot1
[2325]	578
	579
	580	ENDIF
[38]	581
[2325]	582	c----------------------------------------------------------------------
[2844]	583	c Loop on sites:
[2325]	584	c----------------------------------------------------------------------
[38]	585
[2325]	586	c----------------------------------------------------------------------
[2824]	587	c interapolate using temperature in the 7th layer, of surface pressure
[2325]	588	c----------------------------------------------------------------------
[38]	589
[2325]	590	IF(.NOT.firstcal) THEN
	591
[2844]	592	DO iv=1,3
[2325]	593
[38]	594	zp1=0.
	595	zp2=0.
	596	zp2_sm=0.
	597	zt=0.
[2325]	598
[38]	599	DO jj=0,1
[2325]	600
[2844]	601	j=jsite(iv)+jj
[2325]	602
[38]	603	DO ii=0,1
[2325]	604
[2844]	605	i=isite(iv)+ii
[38]	606	zt=zt+zw(ii,jj,iv)*t7(i,j)
	607	zp1=zp1+zw(ii,jj,iv)*log(ps(i,j)) ! interpolate in log(P)
[2824]	608	WRITE (,) 'ps around iv',ps(i,j),iv
[2325]	609
[38]	610	ENDDO
	611	ENDDO
[2325]	612
[38]	613	zp1=exp(zp1) ! because of the bilinear interpolation in log(P)
[2325]	614	WRITE (,) 'constR ',constR
	615	WRITE (,) 'zt ',zt
	616	gh=constR*zt
	617
	618	c----------------------------------------------------------------------
[2824]	619	c surface pressure extrapolated using temp. from 7th atmospheric layer
[2325]	620	c----------------------------------------------------------------------
	621
[38]	622	if (gh.eq.0) then ! if we don't have temperature values
	623	! assume a scale height of 10km
[2844]	624	zp2=zp1exp(-(phisite(iv)-phisim(iv))/(3.731.e4))
[38]	625	else
[2844]	626	zp2=zp1*exp(-(phisite(iv)-phisim(iv))/gh)
[38]	627	endif
[2325]	628
[2844]	629	WRITE (,) 'iv,pstot,zp2, zp1, phisite(iv),phisim(iv),gh'
	630	WRITE (,) iv,pstot*airtot1,zp2,zp1,phisite(iv),phisim(iv),gh
[2824]	631	WRITE(,) "------"
[2325]	632
	633
[2844]	634	c ------Write 3 files (for Vl1, VL2, Insight) --------
[2824]	635	c Sol or ls or both
[2844]	636	c Ps site VLi (i=1,2) or Inisght at true (interpolated) altitude (Pa) (zp2)
	637	c Ps site VLi (i=1,2) or Insight at GCM altitude (Pa) (zp1)
[2325]	638
	639	IF(Time_unit == 1) THEN
	640	WRITE(ifile(iv)+10,'(3e15.5)') dayw,zp2,zp1
	641	ELSEIF (Time_unit == 2) THEN
	642	WRITE(ifile(iv)+10,'(3e15.5)') dayw_ls,zp2,zp1
	643	ELSE
	644	WRITE(ifile(iv)+10,'(4e15.5)') dayw,dayw_ls,zp2,zp1
[2844]	645	ENDIF
	646
	647
[38]	648	ENDDO
[2325]	649
[38]	650	ENDIF
[2844]	651
[2325]	652
[38]	653	firstcal=.false.
	654
[2325]	655
[38]	656	c======================================================================
[2824]	657	c End of loop of variables in the diagfi file
[38]	658	c======================================================================
[2325]	659
[2824]	660	if (sol.ge.unanj-1.e-5) then
	661	! end of year reached (with some roundoff margin)
	662	! increment iyear
	663	iyear=iyear+1
	664	endif
	665
[38]	666	ENDDO
	667
	668	ierr= NF_CLOSE(nid)
	669
[2824]	670	PRINT*,'End of file',nomfich
	671	print*,'Entrer new file name (without trailing .nc)',
	672	&" or return to end"
[38]	673	READ(5,'(a)',err=9999) nomfich
	674
[2325]	675
[2844]	676
[38]	677	c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[2824]	678	c End of loop on the diagfi files
[38]	679	c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	680
	681	ENDDO
	682
[2844]	683	PRINT,'altitude of VL1',.001phis(isite(1),jsite(1))/g
	684	PRINT,'altitude of VL2',.001phis(isite(2),jsite(2))/g
	685	PRINT,'altitude of Ins',.001phis(isite(3),jsite(3))/g
	686
	687	DO iv=1,3
	688	PRINT*,'Site',iv,' i=',isite(iv),'j =',jsite(iv)
[38]	689	WRITE(6,7777)
[2844]	690	s (rlonv(i)*180./pi,i=isite(iv)-1,isite(iv)+2)
[38]	691	print*
[2844]	692	DO j=jsite(iv)-1,jsite(iv)+2
[38]	693	WRITE(6,'(f8.1,10x,5f7.1)')
[2844]	694	s rlatu(j)180./pi,(phis(i,j)/(g1000.),i=isite(iv)-1,
	695	& isite(iv)+2)
[38]	696	ENDDO
	697	print*
	698	print*,'zw'
	699	write(6,'(2(2f10.4/))') ((zw(ii,jj,iv),ii=0,1),jj=0,1)
[2824]	700	print*,'interpolated altitude (km) ',phisim(iv)/1000./g
[38]	701	ENDDO
	702	PRINT*,'R=',r
[2824]	703	9999 PRINT*,'End '
[38]	704
	705	7777 FORMAT ('latitude/longitude',4f7.1)
[2325]	706
[2844]	707	c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	708	c Diurnal Average
	709	c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[2325]	710
[2844]	711	write(,) 'ici'
	712	DO i=1,kyear+1
	713	WRITE(chr2(1:1),'(i1)') i
	714	IF(i.GE.9) WRITE(chr2,'(i2)') i
	715	DO iv=1,3
	716	if (iv==1) then
	717
	718	open(ifile(iv), file = 'ps_VL1_year'//trim(trim(chr2)))
	719	open(ifile(iv)+20, file = 'ps_VL1_year'//trim(chr2)//'_diurnal')
	720	IF(Time_unit == 1) THEN
	721	write(ifile(iv)+20,'(a)') '# Sol , PS at VL1 at true
	722	& (interpolated) altitude (diurnal mean) , PS at VL1 at
	723	& GCM altitude (diurnal mean)'
	724	ELSEIF (Time_unit == 2) THEN
	725	write(ifile(iv)+20,'(a)') '# Ls (deg) , PS at VL1 at
	726	& true (interpolated) altitude (diurnal mean) , PS at VL1
	727	& at GCM altitude (diurnal mean)'
	728	ELSE
	729	write(ifile(iv)+20,'(a)') '# Sol , Ls (deg) , PS at VL1
	730	& at true (interpolated) altitude (diurnal mean) , PS at VL1
	731	& at GCM altitude (diurnal mean)'
	732	ENDIF
	733
	734	elseif (iv==2) then
	735
	736	open(ifile(iv), file = 'ps_VL2_year'//trim(chr2))
	737	open(ifile(iv)+20, file = 'ps_VL2_year'//trim(chr2)//'_diurnal')
	738	IF(Time_unit == 1) THEN
	739	write(ifile(iv)+20,'(a)') '# Sol , PS at VL2 at true
	740	& (interpolated) altitude (diurnal mean) , PS at VL1 at
	741	& GCM altitude (diurnal mean)'
	742	ELSEIF (Time_unit == 2) THEN
	743	write(ifile(iv)+20,'(a)') '# Ls (deg) , PS at VL2
	744	& at true (interpolated) altitude (diurnal mean) , PS at VL1
	745	& at GCM altitude (diurnal mean)'
	746	ELSE
	747	write(ifile(iv)+20,'(a)') '# Sol , Ls (deg) , PS at VL2
	748	& at true (interpolated) altitude (diurnal mean) , PS at VL1
	749	& at GCM altitude (diurnal mean)'
	750	ENDIF
	751
	752	else
	753	open(ifile(iv), file = 'ps_INS_year'//trim(chr2))
	754	open(ifile(iv)+20, file = 'ps_INS_year'//trim(chr2)//'_diurnal')
	755	IF(Time_unit == 1) THEN
	756	write(ifile(iv)+20,'(a)') '# Sol , PS at Insight at true
	757	& (interpolated) altitude (diurnal mean) , PS at VL1
	758	& at GCM altitude (diurnal mean)'
	759	ELSEIF (Time_unit == 2) THEN
	760	write(ifile(iv)+20,'(a)') '# Ls (deg) , PS at Insight at true
	761	& (interpolated) altitude (diurnal mean) , PS at VL1
	762	& at GCM altitude (diurnal mean)'
	763	ELSE
	764	write(ifile(iv)+20,'(a)') '# Sol , Ls (deg) , PS at Insight
	765	& at true (interpolated) altitude (diurnal mean) , PS at VL1
	766	& at GCM altitude (diurnal mean)'
	767	ENDIF
	768	endif
	769
	770	READ(ifile(iv),*)
	771	IF(Time_unit == 1) THEN
	772	READ(ifile(iv),*) dayw,zp2,zp1
	773	ELSEIF (Time_unit == 2) THEN
	774	READ(ifile(iv),*) dayw_ls,zp2,zp1
	775	dayw = ls2sol(dayw_ls)
	776	ELSE
	777	READ(ifile(iv),*) dayw,dayw_ls,zp2,zp1
	778	ENDIF
	779
	780	di=floor(dayw)
	781	di_prev = floor(dayw)
	782
	783	DO k=1,nbmax
	784	ps_gcm_diurnal = 0.
	785	ps_diurnal = 0.
	786	compt_diurn = 0
	787
	788	DO WHILE (di==di_prev)
	789
	790	ps_gcm_diurnal = ps_gcm_diurnal + zp1
	791	ps_diurnal = ps_diurnal + zp2
	792	compt_diurn = compt_diurn + 1
	793	IF(Time_unit == 1) THEN
	794	READ(ifile(iv),*,end=777) dayw,zp2,zp1
	795	ELSEIF (Time_unit == 2) THEN
	796	READ(ifile(iv),*,end=777) dayw_ls,zp2,zp1
	797	dayw=ls2sol(dayw_ls)
	798	ELSE
	799	READ(ifile(iv),*,end=777) dayw,dayw_ls,zp2,zp1
	800	ENDIF
	801	di=floor(dayw)
	802
	803	ENDDO
	804
	805	ps_gcm_diurnal = ps_gcm_diurnal/compt_diurn
	806	ps_diurnal = ps_diurnal/compt_diurn
	807
	808	IF(Time_unit == 1) THEN
	809
	810	write(ifile(iv)+20,'(i4,2e15.5)') di_prev, ps_diurnal,
	811	& ps_gcm_diurnal
	812
	813	ELSEIF (Time_unit == 2) THEN
	814
	815	write(ifile(iv)+20,'(3e15.5)') dayw_ls, ps_diurnal,
	816	& ps_gcm_diurnal
	817
	818	ELSE
	819
	820	write(ifile(iv)+20,'(i4,3e15.5)') di_prev, dayw_ls,
	821	& ps_diurnal, ps_gcm_diurnal
	822
	823	ENDIF
	824
	825
	826	di_prev = di
	827
	828	ENDDO
	829	close(ifile(iv)+20)
	830	close(ifile(iv))
	831	777 ENDDO
	832	ENDDO
[2325]	833
[2844]	834
	835
	836	c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	837	c Harmonics
	838	c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	839
	840
	841	DO i=1,kyear+1
	842	WRITE(chr2(1:1),'(i1)') i
	843	IF(i.GE.9) WRITE(chr2,'(i2)') i
	844	DO iv=1,3
	845
	846	if (iv==1) then
	847	open(ifile(iv), file = 'ps_VL1_year'//trim(chr2)//'_diurnal')
	848	open(ifile(iv)+20, file = 'ps_VL1_year'//trim(chr2)//'_harmonics')
	849
	850	IF(Time_unit == 1) THEN
	851	write(ifile(iv)+20,'(a)') '# Sol , PS at VL1 at true
	852	& (interpolated) altitude (harmonics fit) , PS at VL1 at
	853	& GCM altitude (harmonics fit)'
	854	ELSEIF (Time_unit == 2) THEN
	855	write(ifile(iv)+20,'(a)') '# Ls (deg) , PS at VL1 at
	856	& true (interpolated) altitude (harmonics fit) , PS at VL1
	857	& at GCM altitude (harmonics fit)'
	858	ELSE
	859	write(ifile(iv)+20,'(a)') '# Sol , Ls (deg) , PS at VL1
	860	& at true (interpolated) altitude (harmonics fit) , PS at VL1
	861	& at GCM altitude (harmonics fit)'
	862	ENDIF
	863
	864	elseif (iv==2) then
	865	open(ifile(iv), file = 'ps_VL2_year'//trim(chr2)//'_diurnal')
	866	open(ifile(iv)+20, file = 'ps_VL2_year'//trim(chr2)//'_harmonics')
	867
	868	IF(Time_unit == 1) THEN
	869	write(ifile(iv)+20,'(a)') '# Sol , PS at VL2 at true
	870	& (interpolated) altitude (harmonics fit) , PS at VL2 at
	871	& GCM altitude (harmonics fit)'
	872	ELSEIF (Time_unit == 2) THEN
	873	write(ifile(iv)+20,'(a)') '# Ls (deg) , PS at VL2 at
	874	& true (interpolated) altitude (harmonics fit) , PS at VL2
	875	& at GCM altitude (harmonics fit)'
	876	ELSE
	877	write(ifile(iv)+20,'(a)') '# Sol , Ls (deg) , PS at VL2
	878	& at true (interpolated) altitude (harmonics fit) , PS at VL2
	879	& at GCM altitude (harmonics fit)'
	880	ENDIF
	881
	882	else
	883	open(ifile(iv), file = 'ps_INS_year'//trim(chr2)//'_diurnal')
	884	open(ifile(iv)+20, file = 'ps_INS_year'//trim(chr2)//'_harmonics')
	885
	886	IF(Time_unit == 1) THEN
	887	write(ifile(iv)+20,'(a)') '# Sol , PS at Insight at true
	888	& (interpolated) altitude (harmonics fit) , PS at Insight at
	889	& GCM altitude (harmonics fit)'
	890	ELSEIF (Time_unit == 2) THEN
	891	write(ifile(iv)+20,'(a)') '# Ls (deg) , PS at Insight at
	892	& true (interpolated) altitude (harmonics fit) , PS at Insight
	893	& at GCM altitude (harmonics fit)'
	894	ELSE
	895	write(ifile(iv)+20,'(a)') '# Sol , Ls (deg) , PS at Insight
	896	& at true (interpolated) altitude (harmonics fit) , PS at Insight
	897	& at GCM altitude (harmonics fit)'
	898	ENDIF
	899
	900	endif
	901
	902	READ(ifile(iv),*)
	903
	904	DO k = 1,nbmax
	905	IF (Time_unit == 1) THEN
	906	READ(ifile(iv),*,end=99) di_prev, ps_diurnal_tab(k),
	907	& ps_gcm_diurnal_tab(k)
	908	call sol2ls(real(di_prev),ls_tab(k))
	909	ELSEIF (Time_unit == 2) THEN
	910	READ(ifile(iv),*,end=99) ls_tab(k), ps_diurnal_tab(k),
	911	& ps_gcm_diurnal_tab(k)
	912	ELSE
	913	READ(ifile(iv),*,end=99) di_prev, ls_tab(k),
	914	& ps_diurnal_tab(k), ps_gcm_diurnal_tab(k)
	915	ENDIF
	916
	917	ENDDO
	918
	919
	920	99 ndata=k-1
	921
	922	if(ls_tab(ndata).gt.350.) then
	923
	924	do k = 0,n_harmo
	925
	926	call DiscreetFourierHn(ndata,ls_tab,ps_diurnal_tab,k,a,b)
	927	if(modulo(k,2)==0) then
	928	a_tab(k+1)= a
	929	b_tab(k+1)= b
	930	else
	931	a_tab(k+1)= -a
	932	b_tab(k+1)= -b
	933	endif
	934
	935	call DiscreetFourierHn(ndata,ls_tab,ps_gcm_diurnal_tab,k,a,b)
	936	if(modulo(k,2)==0) then
	937	a_tab_gcm(k+1)= a
	938	b_tab_gcm(k+1)= b
	939	else
	940	a_tab_gcm(k+1)= -a
	941	b_tab_gcm(k+1)= -b
	942	endif
	943
	944	enddo
	945
	946	write(ifile(iv)+20,'(a)') 'Fourrier coefficients ak/bk
	947	&(interpolated altitude)'
	948	write(ifile(iv)+20,'(a,7f)') '#',(a_tab(k),k=1,n_harmo+1)
	949	write(ifile(iv)+20,'(a,7f)') '#',(b_tab(k),k=1,n_harmo+1)
	950	write(ifile(iv)+20,'(a)') 'Fourrier coefficients ak/bk
	951	&(GCM altitude)'
	952	write(ifile(iv)+20,'(a,7f)') '#',(a_tab_gcm(k),k=1,n_harmo+1)
	953	write(ifile(iv)+20,'(a,7f)') '#',(b_tab_gcm(k),k=1,n_harmo+1)
	954
	955
	956	do l = 1,669
	957
	958	call sol2ls(real(l),ls_harm)
	959	ps_harm = a_tab(1)
	960	ps_gcm_harm = a_tab_gcm(1)
	961
	962	do k = 1,n_harmo
	963
	964	ps_harm = ps_harm + a_tab(k+1)*
	965	&cos((pi/180.)k(ls_harm))
	966	&+ b_tab(k+1)sin((pi/180.)k*(ls_harm))
	967
	968	ps_gcm_harm = ps_gcm_harm + a_tab_gcm(k+1)*
	969	&cos((pi/180.)k(ls_harm))
	970	&+ b_tab_gcm(k+1)sin((pi/180.)k*(ls_harm))
	971
	972	enddo
	973
	974
	975	IF(Time_unit == 1) THEN
	976
	977	write(ifile(iv)+20,'(i4,2e15.5)') l, ps_harm,
	978	& ps_gcm_harm
	979
	980	ELSEIF (Time_unit == 2) THEN
	981
	982	write(ifile(iv)+20,'(3e15.5)') ls_harm, ps_harm,
	983	& ps_gcm_harm
	984
	985	ELSE
	986
	987	write(ifile(iv)+20,'(i4,3e15.5)') l,ls_harm, ps_harm,
	988	& ps_gcm_harm
	989
	990	ENDIF
	991
	992	enddo
	993	close(ifile(iv))
	994	close(ifile(iv)+20)
	995
	996	else
	997	write(ifile(iv)+20,'(a)') 'not computed because
	998	& year is not complete'
	999	endif ! (ls.gt.350)
	1000
	1001	ENDDO
	1002	ENDDO
	1003
	1004
	1005
[38]	1006	END
[2325]	1007
	1008	subroutine sol2ls(sol,Ls)
	1009	!==============================================================================
	1010	! Purpose:
	1011	! Convert a date/time, given in sol (martian day),
	1012	! into solar longitude date/time, in Ls (in degrees),
	1013	! where sol=0 is (by definition) the northern hemisphere
	1014	! spring equinox (where Ls=0).
	1015	!==============================================================================
	1016	! Notes:
	1017	! Even though "Ls" is cyclic, if "sol" is greater than N (martian) year,
	1018	! "Ls" will be increased by N*360
	1019	! Won't work as expected if sol is negative (then again,
	1020	! why would that ever happen?)
	1021	!==============================================================================
	1022
	1023	implicit none
	1024
	1025	!==============================================================================
	1026	! Arguments:
	1027	!==============================================================================
	1028	real,intent(in) :: sol
	1029	real,intent(out) :: Ls
	1030
	1031	!==============================================================================
	1032	! Local variables:
	1033	!==============================================================================
	1034	real year_day,peri_day,timeperi,e_elips,twopi,degrad
	1035	data year_day /669./ ! # of sols in a martian year
	1036	data peri_day /485.0/
	1037	data timeperi /1.9082314/
	1038	data e_elips /0.093358/
	1039	data twopi /6.2831853/ ! 2.*pi
	1040	data degrad /57.2957795/ ! pi/180
	1041
	1042	real zanom,xref,zx0,zdx,zteta,zz
	1043
	1044	integer count_years
	1045	integer iter
	1046
	1047	!==============================================================================
	1048	! 1. Compute Ls
	1049	!==============================================================================
	1050
	1051	zz=(sol-peri_day)/year_day
	1052	zanom=twopi*(zz-nint(zz))
	1053	xref=abs(zanom)
	1054
	1055	! The equation zx0 - e * sin (zx0) = xref, solved by Newton
	1056	zx0=xref+e_elips*sin(xref)
	1057	do iter=1,20 ! typically, 2 or 3 iterations are enough
	1058	zdx=-(zx0-e_elipssin(zx0)-xref)/(1.-e_elipscos(zx0))
	1059	zx0=zx0+zdx
	1060	if(abs(zdx).le.(1.e-7)) then
	1061	! write(,)'iter:',iter,' \|zdx\|:',abs(zdx)
	1062	exit
	1063	endif
	1064	enddo
	1065
	1066	if(zanom.lt.0.) zx0=-zx0
	1067
	1068	zteta=2.atan(sqrt((1.+e_elips)/(1.-e_elips))tan(zx0/2.))
	1069	Ls=zteta-timeperi
	1070
	1071	if(Ls.lt.0.) then
	1072	Ls=Ls+twopi
	1073	else
	1074	if(Ls.gt.twopi) then
	1075	Ls=Ls-twopi
	1076	endif
	1077	endif
	1078
	1079	Ls=degrad*Ls
	1080	! Ls is now in degrees
	1081
	1082	!==============================================================================
	1083	! 1. Account for (eventual) years included in input date/time sol
	1084	!==============================================================================
	1085
[2824]	1086	count_years=0 ! initialize
[2325]	1087	zz=sol ! use "zz" to store (and work on) the value of sol
	1088	do while (zz.ge.year_day)
	1089	count_years=count_years+1
	1090	zz=zz-year_day
	1091	enddo
	1092
	1093	! Add 360 degrees to Ls for every year
	1094	if (count_years.ne.0) then
	1095	Ls=Ls+360.*count_years
	1096	endif
	1097
	1098
	1099	end subroutine sol2ls
[2844]	1100
	1101
	1102	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1103	real function ls2sol(ls)
	1104
	1105	! Returns solar longitude, Ls (in deg.), from day number (in sol),
	1106	! where sol=0=Ls=0 at the northern hemisphere spring equinox
	1107
	1108	implicit none
	1109
	1110	! Arguments:
	1111	real, intent(in) :: ls
	1112
	1113	! Local:
	1114	double precision xref,zx0,zteta,zz
	1115	! xref: mean anomaly, zteta: true anomaly, zx0: eccentric anomaly
	1116	double precision year_day
	1117	double precision peri_day,timeperi,e_elips
	1118	double precision pi,degrad
	1119	parameter (year_day=668.6d0) ! number of sols in a amartian year
	1120	! data peri_day /485.0/
	1121	parameter (peri_day=485.35d0) ! date (in sols) of perihelion
	1122	! timeperi: 2pi( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99
	1123	parameter (timeperi=1.90258341759902d0)
	1124	parameter (e_elips=0.0934d0) ! eccentricity of orbit
	1125	parameter (pi=3.14159265358979d0)
	1126	parameter (degrad=57.2957795130823d0)
	1127
	1128	if (abs(ls).lt.1.0e-5) then
	1129	if (ls.ge.0.0) then
	1130	ls2sol = 0.0
	1131	else
	1132	ls2sol = real(year_day)
	1133	end if
	1134	return
	1135	end if
	1136
	1137	zteta = ls/degrad + timeperi
	1138	zx0 = 2.0datan(dtan(0.5zteta)/dsqrt((1.+e_elips)/(1.-e_elips)))
	1139	xref = zx0-e_elips*dsin(zx0)
	1140	zz = xref/(2.*pi)
	1141	ls2sol = real(zz*year_day + peri_day)
	1142	if (ls2sol.lt.0.0) ls2sol = ls2sol + real(year_day)
	1143	if (ls2sol.ge.year_day) ls2sol = ls2sol - real(year_day)
	1144
	1145	return
	1146	end
	1147
	1148	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1149
	1150
	1151	!****************************************************************
	1152	!* Calculate the Fourier harmonic #n of a periodic discreet *
	1153	!* function F(x) defined by ndata points. *
	1154	!* ------------------------------------------------------------ *
	1155	!* Inputs: *
	1156	!* ndata: number of points of discreet function. *
	1157	!* X : pointer to table storing xi abscissas. *
	1158	!* Y : pointer to table storing yi ordinates. *
	1159	!* *
	1160	!* Outputs: *
	1161	!* a : coefficient an of the Fourier series. *
	1162	!* b: : coefficient bn of the Fourier series. *
	1163	!* ------------------------------------------------------------ *
	1164	!* Reference: "Mathematiques en Turbo-Pascal By Marc Ducamp and *
	1165	!* Alain Reverchon, 1. Analyse, Editions Eyrolles, *
	1166	!* Paris, 1991" [BIBLI 03]. *
	1167	!* *
	1168	!* F90 Version By J-P Moreau, Paris. *
	1169	!* (www.jpmoreau.fr) *
	1170	!****************************************************************
	1171	! Note: The Fourier series of a periodic discreet function F(x)
	1172	! can be written under the form:
	1173	! n=inf.
	1174	! F(x) = a0 + Sum ( an cos(2 n pi/T x) + bn sin(2 n pi/T x)
	1175	! n=1
	1176	! ----------------------------------------------------------------
	1177
	1178	Subroutine DiscreetFourierHn(ndata, X, Y, n, a, b)
	1179	real X(1:ndata), Y(1:ndata), a, b
	1180	integer ndata,n,i
	1181	real xi,T,om,wa,wb,wc,wd,wg,wh,wi,wl,wm,wn,wp
	1182	real PI
	1183	PI=4.d0*datan(1.d0)
	1184	T=X(ndata)-X(1); xi=(X(ndata)+X(1))/2.d0
	1185	om = 2PIn/T; a=0.d0; b=0.d0
	1186	do i=1, ndata-1
	1187	wa=X(i); wb=X(i+1)
	1188	wc=Y(i); wd=Y(i+1)
	1189	if (wa.ne.wb) then
	1190	wg = (wd-wc)/(wb-wa)
	1191	wh = om(wa-xi); wi=om(wb-xi)
	1192	if (n==0) then
	1193	a = a + (wb-wa)(wc+wg/2.d0(wb-wa))
	1194	else
	1195	wl = cos(wh); wm = sin(wh)
	1196	wn = cos(wi); wp = sin(wi)
	1197	a = a + wg/om(wn-wl) + wdwp - wc*wm
	1198	b = b + wg/om(wp-wm) - wdwn + wc*wl
	1199	end if
	1200	end if
	1201	end do
	1202	a = a/T; b = b/T
	1203	if (n.ne.0) then
	1204	a = a2.d0/om; b = b2.d0/om
	1205	end if
	1206	return
	1207	End

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

Download in other formats:

Original Format