Context Navigation

fit_Iceinertia_MONSicedepth.F @ 2970

Last change on this file since 2970 was 2567, checked in by emillour, 3 years ago
Mars GCM utilities: Minor fixes to run with picky gfortran 10.3.0 which requires one element arrays (rather than scalars) when calling NetCDF routines, andf that stop statements should not be followed by strings. While at it replaced tabs with spaces. EM
File size: 34.3 KB

Rev	Line
[2338]	1	PROGRAM SUPERFIT
	2
	3	c ------------------------------------------------------------------
	4	c #### Special version with ice depth and ice inertia #####
	5	c Program used to compute best-fit cap albedo, ICE DEPTH and total CO2
	6	c inventory
	7	c Make use of 2 annual GCM run performed with 2 different cap albedos.
	8	c Based on Hourdin et al. (JGR, 1995)
	9	c Modified to account for non-linear sensitivity of polar cap mass
	10	c To cap albedo (FF,2000)
	11	c Modified/cleaned-up to use 4 input files and minimise cap albedos,
	12	c MONS-derived ice depth and total pressure (EM,2008)
	13	c Modified to minimise wrt ice thermal inertia coefficients (EM, 2009)
	14	c ------------------------------------------------------------------
	15
	16	implicit none
	17
	18	integer ngcm,nvl1,nsol,time_unit
	19
	20	c *** THE FOLLOWING LINES MUST BE ADAPTED FOR EACH CASES : ************
	21	! parameter(ngcm=10704) ! size of the 1 year raw GCM data
	22	! parameter(ngcm=2672) ! "4 & 6" files
	23	parameter(ngcm=8028) ! "3 & 5" files
	24	parameter(nvl1=669) ! size of the Vl1 Ps Observation
	25	c *************************************************************
	26
	27
	28	c size the smoothed data (1/sol) which are used for the simulation
	29	parameter(nsol=669) ! size the smoothed data (1/sol)
	30	c parameter(nsol=445) ! size the smoothed data (1/sol)
	31
	32	c OBSERVED VL1 pressure
	33	c ~~~~~~~~~~~~~~~~~~~~~
	34	real pvl_obs(nvl1), solvl(nvl1), lsvl(nvl1)
	35
	36	c Raw Data from 4 simulations with 2 different albedo/iced :
	37	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	38	c runs are numbered 1-4 such than runs 1&2 are at same ice depth
	39	c and runs 3&4 are at same ice depth
	40	c and runs 1&3 are at same ice inertia
	41	c and runs 2&4 are at same ice inertia
	42	real iceigcmN(4) , iceigcmS(4)
	43	real icedgcm(4)
	44	integer refrun ! run number of the 'reference' run
	45	! reference values for functions (e.g. values of run #1)
	46	real iceigcm_refN, icedgcm_ref
	47	real iceigcm_refS
	48
	49	real ptot(4)
	50	real patm(ngcm,4), pcapn(ngcm,4),pcaps(ngcm,4)
	51	real pvl_gcm(ngcm,4) ! Simulated VL1 pressure
	52	real solgcm(ngcm) ! time in runs in sol
	53	real lsgcm(ngcm) ! time in runs in ls
	54	integer year_xvik ! year of xvik files to use for fit
	55	integer plot_test ! =1 if output file is to plot COST(DN,DS), =2 for COST(IN,IS)s
	56
	57	c Smoothed Data from the 2 simulations with 2 different albedo :
	58	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	59	real box
	60	integer n ,i
	61	real sol(nsol) ! time in runs
	62	real ls(nsol)
	63	real solconv
	64	real pvl_sm(nsol,4) ! Simulated smooth VL1 pressure
	65	real patm_sm(nsol,4), pcapn_sm(nsol,4),pcaps_sm(nsol,4)
	66	real alphavl1(nsol) ! pvl/patm ratio
	67	c Mass cap sensitivity to thermal inertia (derivative)
	68	real dpdicein(nsol),dpdiceis(nsol)
	69	real dpdicein_fltr(nsol),dpdiceis_fltr(nsol)
	70	c Mass cap sensitivity to ice inertia (derivative)
	71	real dpden(nsol),dpdes(nsol)
	72	real dpden_fltr(nsol),dpdes_fltr(nsol)
	73
	74	real iceis, icein, ptry ,pvl1
	75	real iceds, icedn
	76	! ice inertia range to explore, and maximum allowed difference between N&S ice inertias:
	77	real iceimin,iceimax,maxiceidiff
	78	real deltaicei ! step in ice thermal inertia
	79	! ice depth range to explore, and maximum allowed difference between N&S ice depths:
	80	real icedmin,icedmax,maxiceddiff
	81	real deltaiced ! step in ice depth coefficient
	82	! total pressure range to explore, and pressure step:
	83	real pmin,pmax,deltap
	84	real cost,cost4plot
	85	real iceis4plot,ps4plot,icein4plot !,iceds4plot
	86	real iceds4plot,icedn4plot
	87	real costmin, pfit, iceinfit, iceisfit
	88	real pcapn_new, pcaps_new
	89	real icednfit, icedsfit
	90	logical fit
	91
	92	real fonc, fonc2n, fonc2s
	93
	94	c Inputs files (from xvik) variables :
	95	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	96	character(len=33) filename1, filename2
	97	character(len=99) dset_DminImin,dset_DmaxImin
	98	character(len=99) dset_DminImax,dset_DmaxImax
	99	integer ie
	100
	101
	102	c----------------------------------------------------------------------
	103	c Initialisation :
	104	c ~~~~~~~~~~~~~~
	105	c Implicit function : behavior of Pcap = fonc(alb)
	106	c (e.g. fonc(alb) = alb if linear)
	107	!! fonc(albn) = exp(3.46*albn) ! for albedo, both north & south
	108	fonc(icein) = icein
	109	c fonc2n(icedn) = log(icedn +10)
	110	fonc2n(icedn) = icedn ! for northern ice depth coefficient
	111	c fonc2s(iceds) = log(iceds +1)
	112	fonc2s(iceds) = iceds ! for southern ice depth coefficient
	113
	114	c *** THE FOLLOWING LINES MUST BE ADAPTED FOR EACH CASES : ************
	115
	116	write(,) 'Program written for xvik outputs files ''xpsol'' ',
	117	&'and ''xprestot'' from xvik.F'
	118	write(,) 'Enter extreme parameters values used in xvik.F'
	119	write(,) 'Minimum ice depth coefficient'
	120	read(,) icedmin
	121	write(,) icedmin
	122	write(,) 'Maximum ice depth coefficient'
	123	read(,) icedmax
	124	write(,) icedmax
	125	write(,) 'Minimum ice thermal inertia '
	126	read(,) iceimin
	127	write(,) iceimin
	128	write(,) 'Maximum ice thermal inertia '
	129	read(,) iceimax
	130	write(,) iceimax
	131
	132
	133	maxiceidiff = 3000 ! maximum allowed difference between N and S best fit ice thermal inertia
[2567]	134	deltaicei= 20 ! step in ice thermal inertia
[2338]	135	maxiceddiff = 30.e-4 ! maximum allowed difference between N and S best fit ice depth
	136	deltaiced= 0.5e-4 ! step in ice depth coefficient
	137	pmin = 690 ! minimum allowed best fit CO2 total pressure (Pa)
	138	pmax = 720 ! maximum allowed best fit CO2 total pressure (Pa)
	139	deltap=1.0 ! step in pressure
	140
	141	c File to be read :
	142	c ~~~~~~~~~~~~~~~
	143	c runs are numbered 1-4 such than runs 1&2 are at same ice depth
	144	c and runs 3&4 are at same ice depth
	145	c and runs 1&3 are at same ice inertia
	146	c and runs 2&4 are at same ice inertia
	147	! First run:
	148	iceigcmN(1)= iceimin ! Northern Cap ice thermal inertia
	149	iceigcmS(1)= iceimin ! Southern Cap ice thermal inertia
	150	icedgcm(1) = icedmin ! Ice depth coefficient
	151	! Second run:
	152	iceigcmN(2) = iceimax ! Northern Cap ice thermal inertia
	153	iceigcmS(2) = iceimax ! Southern Cap ice thermal inertia
	154	icedgcm(2) = icedmin ! Ice depth coefficient
	155	! Third run:
	156	iceigcmN(3) = iceimin ! Northern Cap ice thermal inertia
	157	iceigcmS(3) = iceimin ! Southern Cap ice thermal inertia
	158	icedgcm(3) = icedmax ! Ice depth coefficient
	159	! Fourth run:
	160	iceigcmN(4) = iceimax ! Northern Cap ice thermal inertia
	161	iceigcmS(4) = iceimax ! Southern Cap ice thermal inertia
	162	icedgcm(4) = icedmax ! Ice depth coefficient
	163
	164	! set reference values: those of one of the files (here file 1)
	165	refrun=1
	166	! ice inertia of the reference run used to simulate function
	167	iceigcm_refN=iceigcmN(refrun)
	168	! ice inertia of the reference run used to simulate function
	169	iceigcm_refS=iceigcmS(refrun)
	170	! icedepth coefficient of the reference run used to simulate function
	171	icedgcm_ref=icedgcm(refrun)
	172
	173	write(,) 'Path to xvik outputs with minimum ice depth ',
	174	&'coefficient and minimum ice thermal inertia ?'
	175	read (*,'(a)') dset_DminImin
	176	write(,) dset_DminImin
	177	write(,) 'Path to xvik outputs with minimum ice depth ',
	178	&'coefficient and maximum ice thermal inertia ?'
	179	read (*,'(a)') dset_DminImax
	180	write(,) dset_DminImax
	181	write(,) 'Path to xvik outputs with maximum ice depth ',
	182	&'coefficient and minimum ice thermal inertia ?'
	183	read (*,'(a)') dset_DmaxImin
	184	write(,) dset_DmaxImin
	185	write(,) 'Path to xvik outputs with maximum ice depth ',
	186	&'coefficient and maximum ice thermal inertia ?'
	187	read (*,'(a)') dset_DmaxImax
	188	write(,) dset_DmaxImax
	189
	190	write(,) 'Which year of xvik outputs do you want to use ?'
	191	read (,) year_xvik
	192	write(,) year_xvik
	193
	194	write(,) 'Xvik files in sol (1), ls (2) or both (3)'
	195	read(,) time_unit
	196	write(,) time_unit
	197
	198
	199	write(,) 'COST being the model/obs difference'
	200	write(,) 'IN ans IS being the ice thermal inertia of ',
	201	&'Northern and Southern Cap'
	202	write(,) 'DN ans DS being the ice depth coefficient of ',
	203	&'of Northern and Southern Cap'
	204	write(,) 'Do you want to use output file ',
	205	&'''minimization.txt'' to plot COST(DN,DS) (1) or COST(IN,IS) (2)'
	206	read(,) plot_test
	207	write(,) plot_test
	208
	209	c For four simulation, files : "time (VL1 sol) , Ps Vl1 (Pa)"
	210	c "time (VL1 sol), Patm,PcapN,PcapS (Pa)"
	211
	212
	213
	214	write(filename1,fmt='(a6,i1)') 'xpsol1',year_xvik
	215	write(filename2,fmt='(a8,i1)') 'xprestot',year_xvik
	216
	217
	218	! Dmin Imin files
	219
	220	write(,) 'Opening ', trim(dset_DminImin)//'/'//trim(filename1)
	221	open(21,file=trim(dset_DminImin)//'/'//trim(filename1),iostat=ie)
	222
	223	if (ie.ne.0) then
	224	write(,) 'Error opening file ',trim(dset_DminImin)//'/'
	225	&//trim(filename1)
	226	stop
	227	endif
	228
	229	write(,) 'Opening ', trim(dset_DminImin)//'/'//trim(filename2)
	230	open(11,file=trim(dset_DminImin)//'/'//trim(filename2),iostat=ie)
	231
	232	if (ie.ne.0) then
	233	write(,) 'Error opening file ',trim(dset_DminImin)//'/'
	234	&//trim(filename2)
	235	stop
	236	endif
	237
	238	! Dmin Imax files
	239
	240	write(,) 'Opening ', trim(dset_DminImax)//'/'//trim(filename1)
	241	open(22,file=trim(dset_DminImax)//'/'//trim(filename1),iostat=ie)
	242
	243	if (ie.ne.0) then
	244	write(,) 'Error opening file ',trim(dset_DminImax)//'/'
	245	&//trim(filename1)
	246	stop
	247	endif
	248
	249	write(,) 'Opening ', trim(dset_DminImax)//'/'//trim(filename2)
	250	open(12,file=trim(dset_DminImax)//'/'//trim(filename2),iostat=ie)
	251
	252	if (ie.ne.0) then
	253	write(,) 'Error opening file ',trim(dset_DminImax)//'/'
	254	&//trim(filename2)
	255	stop
	256	endif
	257
	258	! Dmax Imin files
	259
	260	write(,) 'Opening ', trim(dset_DmaxImin)//'/'//trim(filename1)
	261	open(23,file=trim(dset_DmaxImin)//'/'//trim(filename1),iostat=ie)
	262
	263	if (ie.ne.0) then
	264	write(,) 'Error opening file ',trim(dset_DmaxImin)//'/'
	265	&//trim(filename1)
	266	stop
	267	endif
	268
	269	write(,) 'Opening ', trim(dset_DmaxImin)//'/'//trim(filename2)
	270	open(13,file=trim(dset_DmaxImin)//'/'//trim(filename2),iostat=ie)
	271
	272	if (ie.ne.0) then
	273	write(,) 'Error opening file ',trim(dset_DmaxImin)//'/'
	274	&//trim(filename2)
	275	stop
	276	endif
	277
	278
	279	! Dmax Imax files
	280
	281	write(,) 'Opening ', trim(dset_DmaxImax)//'/'//trim(filename1)
	282	open(24,file=trim(dset_DmaxImax)//'/'//trim(filename1),iostat=ie)
	283
	284	if (ie.ne.0) then
	285	write(,) 'Error opening file ',trim(dset_DmaxImax)//'/'
	286	&//trim(filename1)
	287	stop
	288	endif
	289
	290	write(,) 'Opening ', trim(dset_DmaxImax)//'/'//trim(filename2)
	291	open(14,file=trim(dset_DmaxImax)//'/'//trim(filename2),iostat=ie)
	292
	293	if (ie.ne.0) then
	294	write(,) 'Error opening file ',trim(dset_DmaxImax)//'/'
	295	&//trim(filename2)
	296	stop
	297	endif
	298
	299
	300
	301	c *************************************************************
	302	c Observed Smooth Viking Curves
	303	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	304	open(30,file='VL1')
	305	do n=1,nsol
	306	read(30,*) solvl(n), lsvl(n), pvl_obs(n)
	307	end do
	308	close(30)
	309
	310	c Opening output file
	311
	312	!open(33, file = 'minimization.txt')
	313
	314
	315	c reading Simulation results
	316	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~
	317
	318
	319	do n=1,ngcm
	320
	321	c Reading Viking Lander 1 simulated pressure for 4 runs
	322
	323	if (time_unit == 1) then
	324	read(21,*) solgcm(n), pvl_gcm(n,1)
	325	read(22,*) solgcm(n), pvl_gcm(n,2)
	326	read(23,*) solgcm(n), pvl_gcm(n,3)
	327	read(24,*) solgcm(n), pvl_gcm(n,4)
	328
	329	c Reading atmospheric and "cap" total pressure for all runs
	330	read(11,*) solgcm(n), patm(n,1), pcapn(n,1),pcaps(n,1)
	331	read(12,*) solgcm(n), patm(n,2), pcapn(n,2),pcaps(n,2)
	332	read(13,*) solgcm(n), patm(n,3), pcapn(n,3),pcaps(n,3)
	333	read(14,*) solgcm(n), patm(n,4), pcapn(n,4),pcaps(n,4)
[2567]	334
[2338]	335	elseif (time_unit == 2) then
	336	read(21,*) lsgcm(n), pvl_gcm(n,1)
	337	read(22,*) lsgcm(n), pvl_gcm(n,2)
	338	read(23,*) lsgcm(n), pvl_gcm(n,3)
	339	read(24,*) lsgcm(n), pvl_gcm(n,4)
	340
	341	c Reading atmospheric and "cap" total pressure for all runs
	342	read(11,*) lsgcm(n), patm(n,1), pcapn(n,1),pcaps(n,1)
	343	read(12,*) lsgcm(n), patm(n,2), pcapn(n,2),pcaps(n,2)
	344	read(13,*) lsgcm(n), patm(n,3), pcapn(n,3),pcaps(n,3)
	345	read(14,*) lsgcm(n), patm(n,4), pcapn(n,4),pcaps(n,4)
	346	call ls2sol(lsgcm(n),solconv)
	347	solgcm(n) = solconv
	348
	349	elseif (time_unit == 3) then
	350	read(21,*) solgcm(n), lsgcm(n), pvl_gcm(n,1)
	351	read(22,*) solgcm(n), lsgcm(n), pvl_gcm(n,2)
	352	read(23,*) solgcm(n), lsgcm(n), pvl_gcm(n,3)
	353	read(24,*) solgcm(n), lsgcm(n), pvl_gcm(n,4)
	354
	355	c Reading atmospheric and "cap" total pressure for all runs
	356	read(11,*) solgcm(n), lsgcm(n), patm(n,1), pcapn(n,1),pcaps(n,1)
	357	read(12,*) solgcm(n), lsgcm(n), patm(n,2), pcapn(n,2),pcaps(n,2)
	358	read(13,*) solgcm(n), lsgcm(n), patm(n,3), pcapn(n,3),pcaps(n,3)
	359	read(14,*) solgcm(n), lsgcm(n), patm(n,4), pcapn(n,4),pcaps(n,4)
[2567]	360
	361	else
	362	write(,) 'Wrong integer for xvik files format :',
[2338]	363	&' must be 1, 2 or 3'
	364	stop
	365
	366	endif
[2567]	367
[2338]	368	c Checking total CO2 inventory for all runs :
	369	do i=1,4
	370	if(n.eq.1) then
	371	ptot(i) = patm(n,i)+pcapn(n,i)+pcaps(n,i)
	372	write(,) 'For run = ',i,' Ptot= ',ptot(i)
	373	else
	374	if(abs(patm(n,i)+pcapn(n,i)+pcaps(n,i)-ptot(i)).gt.3)then
	375	write(,)'total pressure not constant for run i= ',i
	376	write(,) 'n=',1,' ptot=',ptot(i)
	377	write(,)'n=',n,' ptot=',patm(n,i)+pcapn(n,i)+pcaps(n,i)
	378	end if
	379	end if
	380	end do
	381	end do
	382
[2567]	383	close(11)
	384	close(12)
	385	close(13)
	386	close(14)
	387	close(21)
	388	close(22)
	389	close(23)
	390	close(24)
	391
	392
[2338]	393	c Smoothing simulated GCM Viking 1 pressure curves
	394	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	395	do n=1,nsol
	396	sol(n)=float(n)
	397	end do
	398
	399
	400	c Running average for both runs (with "boxsize" box in days)
	401	box=20.
	402	do i=1,4
	403	call runave(solgcm,pvl_gcm(1,i),ngcm,669.,box,
	404	& sol,pvl_sm(1,i),nsol)
	405	call runave(solgcm,patm(1,i),ngcm,669.,box,
	406	& sol,patm_sm(1,i),nsol)
	407	call runave(solgcm,pcapn(1,i),ngcm,669.,box,
	408	& sol,pcapn_sm(1,i),nsol)
	409	call runave(solgcm,pcaps(1,i),ngcm,669.,box,
	410	& sol,pcaps_sm(1,i),nsol)
	411	end do
	412
	413
	414	c Computing mass cap sensitivity to albedo (derivative) and alphaVl1
	415	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	416	c alphaVL1 = Pvl1/patm at a given time (see Hourdin et al., JGR, 1995)
	417
	418	do n=1,nsol
	419	! Evaluate derivative of pressure wrt Northern inertia, as the
	420	! mean of the 2 values obtained using runs 1&2 and 3&4
	421	dpdicein(n)=(((pcapn_sm(n,1)-pcapn_sm(n,2))/
	422	& (fonc(iceigcmN(1))-fonc(iceigcmN(2))))+
	423	& ((pcapn_sm(n,3)-pcapn_sm(n,4))/
	424	& (fonc(iceigcmN(3))-fonc(iceigcmN(4)))))*(1./2.)
	425	! Evaluate derivative of pressure wrt Southern inertia, as the
	426	! mean of the 2 values obtained using runs 1&2 and 3&4
	427	dpdiceis(n)=(((pcaps_sm(n,1)-pcaps_sm(n,2))/
	428	& (fonc(iceigcmS(1))-fonc(iceigcmS(2))))+
	429	& ((pcaps_sm(n,3)-pcaps_sm(n,4))/
	430	& (fonc(iceigcmS(3))-fonc(iceigcmS(4)))))*(1./2.)
	431	! Evaluate derivative of pressure wrt Northern ice depth coefficient,
	432	! as the mean of the 2 values obtained using runs 1&3 and 2&4
	433	dpden(n)=(((pcapn_sm(n,1)-pcapn_sm(n,3))/
	434	& (fonc2n(icedgcm(1))-fonc2n(icedgcm(3))))+
	435	& ((pcapn_sm(n,2)-pcapn_sm(n,4))/
	436	& (fonc2n(icedgcm(2))-fonc2n(icedgcm(4)))))*(1./2.)
	437	! Evaluate derivative of pressure wrt Southern ice depth coefficient,
	438	! as the mean of the 2 values obtained using runs 1&3 and 2&4
	439	dpdes(n)=(((pcaps_sm(n,1)-pcaps_sm(n,3))/
	440	& (fonc2s(icedgcm(1))-fonc2s(icedgcm(3))))+
	441	& ((pcaps_sm(n,2)-pcaps_sm(n,4))/
	442	& (fonc2s(icedgcm(2))-fonc2s(icedgcm(4)))))*(1./2.)
	443	! Evaluate alphaVL1 coefficient, as the mean of alphaVL1 coefficients
	444	! of all 4 runs:
	445	alphavl1(n) = (1./4.)*(pvl_sm(n,1)/patm_sm(n,1)
	446	& +pvl_sm(n,2)/patm_sm(n,2)+pvl_sm(n,3)/patm_sm(n,3)
	447	& +pvl_sm(n,4)/patm_sm(n,4))
	448	!write(91,*) sol(n), dpden(n), dpdes(n)
	449	c write(,) 'pcapn_sm(n,1),pcapn_sm(n,2)'
	450	c & ,pcapn_sm(n,1),pcapn_sm(n,2)
	451	c write(,)'fonc(albgcmN(1)),fonc(albgcmS(2))',
	452	c & fonc(albgcmN(1)),fonc(albgcmS(2))
	453	c write(,)'albgcmN(1)),albgcmS(2)',
	454	c & albgcmN(1),albgcmS(2)
	455	!write(90,*)pvl_sm(n,1)/patm_sm(n,1),pvl_sm(n,2)/patm_sm(n,2),
	456	c & pvl_sm(n,3)/patm_sm(n,3),pvl_sm(n,4)/patm_sm(n,4)
	457	end do ! of do n=1,nsol
	458
	459	c -------------------------------------------------
	460	c Smooth the derivative like Frederic did ? :
	461	call runave(sol,dpdicein,nsol,669.,60.,
	462	& sol,dpdicein_fltr,nsol)
	463	call runave(sol,dpdiceis,nsol,669.,60.,
	464	& sol,dpdiceis_fltr,nsol)
	465	call runave(sol,dpden,nsol,669.,60.,
	466	& sol,dpden_fltr,nsol)
	467	call runave(sol,dpdes,nsol,669.,60.,
	468	& sol,dpdes_fltr,nsol)
	469	do n=1,nsol
	470	dpdiceis(n)=dpdiceis_fltr(n)
	471	dpdicein(n)=dpdicein_fltr(n)
	472	dpdes(n)=dpdes_fltr(n)
	473	dpden(n)=dpden_fltr(n)
	474	!write(92,*) sol(n), dpden_fltr(n), dpdes_fltr(n)
	475	end do
	476	c -------------------------------------------------
	477
	478
	479	c Compute best fit parameters by minimizing Cost function
	480	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	481	c STUPID ultra-robust minimization method: waisting CPU to save human time...
	482
	483	open(33, file = 'minimization.txt')
	484	c 'minimzation.txt' will contain data to plot COST(DN,DS) or COST(IN,IS) depending on plot_test
	485
	486
	487	c to plot COST(DN,DS) (plot_test=1), loops must be in this order : DN->DS->IN->IS
	488
	489	costmin = 1.e30 ! initialization
	490
	491
	492	if(plot_test.eq.1) then
	493	c to plot COST(DN,DS) (plot_test=1), loops must be in this order : DN->DS->IN->IS
	494	write(,)' DN DS cost IN IS Ps'
[2567]	495	icedn=icedmin ! initialization
[2338]	496	do while (icedn.le.icedmax) ! loop on northern ice depth coefficient
	497	! albn=albmin ! initialization
	498	! do while (albn.le.albmax) ! loop on northern cap albedo
	499	iceds=max(icedn-maxiceddiff,icedmin) ! initialization
	500	do while (iceds.le.min(icedn+maxiceddiff,icedmax))
	501	cost4plot = 1.e30 ! initializationqsub
	502	! iceds=max(icedn-maxiceddiff,icedmin) ! initialization
	503	! do while (iceds.le.min(icedn+maxiceddiff,icedmax))
	504	icein=iceimin ! initialization
	505	do while (icein.le.iceimax) ! loop on northern ice inertia
	506	iceis=max(icein-maxiceidiff,iceimin) ! initialization
	507	do while (iceis.le.min(icein+maxiceidiff,iceimax))
	508	ptry=pmin ! initialization
	509	do while (ptry.le.pmax) ! loop on total pressure
	510	cost =0. !initialization
	511	do n=1,nsol
[2567]	512	!write(,) n
[2338]	513	! Pressure corresponding to Northern cap
	514	pcapn_new=pcapn_sm(n,refrun)+
	515	& (fonc(icein) -fonc(iceigcm_refN)) * dpdicein(n)
	516	& + (fonc2n(icedn) -fonc2n(icedgcm_ref)) * dpden(n)
	517	pcapn_new= max(pcapn_new,0.)
	518	! Pressure corresponding to Southern Cap
	519	pcaps_new=pcaps_sm(n,refrun)+
	520	& (fonc(iceis) -fonc(iceigcm_refS)) * dpdiceis(n)
	521	& + (fonc2s(iceds) -fonc2s(icedgcm_ref)) * dpdes(n)
	522	pcaps_new= max(pcaps_new,0.)
	523	! Pressure at VL1 site
	524	pvl1 = alphavl1(n)*
	525	& (ptry - pcapn_new - pcaps_new)
	526	! cumulative squared differences between predicted
	527	! and observed pressures at VL1 site
	528	cost= cost+ ( pvl_obs(n) - pvl1)**2
[2567]	529
[2338]	530	end do ! of do n=1,nsol
	531	! store parameters which lead to minimum cost
	532	! (i.e. best fit so far)
	533	if(cost.lt.costmin) then
	534	costmin = cost
	535	pfit=ptry
	536	iceinfit=icein
	537	iceisfit=iceis
	538	icednfit=icedn
	539	icedsfit=iceds
	540	end if
[2567]	541
[2338]	542	if(cost.lt.cost4plot) then
	543	c RMS, best pressure and best albedos for these icedsivities value
	544	cost4plot = cost
	545	iceis4plot=iceis
	546	! iceds4plot=iceds
	547	icein4plot=icein
	548	ps4plot=ptry
	549	end if
	550	ptry=ptry+deltap ! increment ptry
	551	end do ! of do while (ptry.le.pmax)
	552	iceis=iceis+deltaicei ! increment iceis
	553	end do ! of do while (iceis.le.min(icein+maxiceidiff,iceimax))
	554	icein=icein+deltaicei ! increment icein
	555	enddo ! of do while (icein.le.iceimax)
	556	! iceds=iceds+deltaiced ! increment iceds
	557	! end do ! of do while (iceds.le.min(icedn+maxiceddiff,icedmax))
	558	! write(*,fmt='(1pe9.2,f5.2,f9.3,1pe9.2,
	559	! & f5.2,f7.2)')
	560	write(*,fmt='(6(1pe10.3))') ! output to screen
	561	& icedn,iceds,sqrt(cost4plot/float(nsol)),
	562	& icein4plot,iceis4plot,ps4plot
	563	write(33,fmt='(6(1pe10.3))') ! output to file
	564	& icedn,iceds,sqrt(cost4plot/float(nsol)),
	565	& icein4plot,iceis4plot,ps4plot
	566	iceds=iceds+deltaiced ! increment iceds
	567	enddo ! of do while (iceds.le.min(icedn+maxiceddiff,icedmax))
	568	write(33,*)' ' ! blank line in output file
	569	! albn=albn+deltaalb ! increment albn
	570	! end do ! of do while (albn.le.albmax)
	571	icedn=icedn+deltaiced
	572	end do ! do while (icedn.le.icedmax)
	573
	574	elseif(plot_test.eq.2) then
	575	c to plot COST(IN,IS) (plot_test=2), loops must be in this order : IN->IS->DN->DS
	576	write(,)' DN DS cost IN IS Ps'
[2567]	577	write(33,*)' DN DS cost IN IS Ps'
	578	icein=iceimin ! initialization
[2338]	579	do while (icein.le.iceimax) ! loop on northern ice depth coefficient
	580	! albn=albmin ! initialization
	581	! do while (albn.le.albmax) ! loop on northern cap albedo
	582	iceis=max(icein-maxiceidiff,iceimin)! initialization
	583	do while (iceis.le.min(icein+maxiceidiff,iceimax))
	584	cost4plot = 1.e30 ! initializationqsub
	585	! iceds=max(icedn-maxiceddiff,icedmin) ! initialization
	586	! do while (iceds.le.min(icedn+maxiceddiff,icedmax))
	587	icedn=icedmin ! initialization
	588	do while (icedn.le.icedmax) ! loop on northern ice inertia
	589	iceds=max(icedn-maxiceddiff,icedmin) ! initialization
	590	do while (iceds.le.min(icedn+maxiceddiff,icedmax))
	591	ptry=pmin ! initialization
	592	do while (ptry.le.pmax) ! loop on total pressure
	593	cost =0. !initialization
	594	do n=1,nsol
[2567]	595	!write(,) n
[2338]	596	! Pressure corresponding to Northern cap
	597	pcapn_new=pcapn_sm(n,refrun)+
	598	& (fonc(icein) -fonc(iceigcm_refN)) * dpdicein(n)
	599	& + (fonc2n(icedn) -fonc2n(icedgcm_ref)) * dpden(n)
	600	pcapn_new= max(pcapn_new,0.)
	601	! Pressure corresponding to Southern Cap
	602	pcaps_new=pcaps_sm(n,refrun)+
	603	& (fonc(iceis) -fonc(iceigcm_refS)) * dpdiceis(n)
	604	& + (fonc2s(iceds) -fonc2s(icedgcm_ref)) * dpdes(n)
	605	pcaps_new= max(pcaps_new,0.)
	606	! Pressure at VL1 site
	607	pvl1 = alphavl1(n)*
	608	& (ptry - pcapn_new - pcaps_new)
	609	! cumulative squared differences between predicted
	610	! and observed pressures at VL1 site
	611	cost= cost+ ( pvl_obs(n) - pvl1)**2
[2567]	612
[2338]	613	end do ! of do n=1,nsol
	614	! store parameters which lead to minimum cost
	615	! (i.e. best fit so far)
	616	if(cost.lt.costmin) then
	617	costmin = cost
	618	pfit=ptry
	619	iceinfit=icein
	620	iceisfit=iceis
	621	icednfit=icedn
	622	icedsfit=iceds
	623	end if
[2567]	624
[2338]	625	if(cost.lt.cost4plot) then
	626	c RMS, best pressure and best albedos for these icedsivities value
	627	cost4plot = cost
	628	iceds4plot=iceds
	629	! iceds4plot=iceds
	630	icedn4plot=icedn
	631	ps4plot=ptry
	632	end if
	633	ptry=ptry+deltap ! increment ptry
	634	end do ! of do while (ptry.le.pmax)
	635	iceds=iceds+deltaiced ! increment iceis
	636	end do ! of do while (iceis.le.min(icein+maxiceidiff,iceimax))
	637	icedn=icedn+deltaiced ! increment icein
	638	enddo ! of do while (icein.le.iceimax)
	639	! iceds=iceds+deltaiced ! increment iceds
	640	! end do ! of do while (iceds.le.min(icedn+maxiceddiff,icedmax))
	641	! write(*,fmt='(1pe9.2,f5.2,f9.3,1pe9.2,
	642	! & f5.2,f7.2)')
	643	write(*,fmt='(6(1pe10.3))') ! output to screen
	644	& icedn4plot,iceds4plot,sqrt(cost4plot/float(nsol)),
	645	& icein,iceis,ps4plot
	646	write(33,fmt='(6(1pe10.3))') ! output to file
	647	& icedn4plot,iceds4plot,sqrt(cost4plot/float(nsol)),
	648	& icein,iceis,ps4plot
	649	iceis=iceis+deltaicei ! increment iceds
	650	enddo ! of do while (iceds.le.min(icedn+maxiceddiff,icedmax))
	651	write(33,*)' ' ! blank line in output file
	652	! albn=albn+deltaalb ! increment albn
	653	! end do ! of do while (albn.le.albmax)
	654	icein=icein+deltaicei
	655	end do ! do while (icedn.le.icedmax)
	656	else
	657	write(,) 'Wrong integer for plot (must be 1 or 2)'
	658	stop
	659
	660	endif
	661
	662	close(33)
	663
	664	write(,) 'Best fit Ptot=', pfit
	665	write(,) 'Best fit In=', iceinfit
	666	write(,) 'Best fit IS=', iceisfit
	667	write(,) 'Best fit Dn=', icednfit
	668	write(,) 'Best fit DS=', icedsfit
	669	write(,) 'RMS difference model/obs=',sqrt(costmin/float(nsol))
	670
	671	c Synthethic VL1 pressure curves for information
	672	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	673	open(41,file='pvlfit')
	674	open(42,file='xprestotfit')
	675	do n=1,nsol
	676	pcapn_new=pcapn_sm(n,refrun)+
	677	& (fonc(iceinfit) -fonc(iceigcm_refN)) * dpdicein(n)
	678	& + (fonc2n(icednfit) -fonc2n(icedgcm_ref)) * dpden(n)
	679	pcapn_new= max(pcapn_new,0.)
	680	pcaps_new=pcaps_sm(n,refrun)+
	681	& (fonc(iceisfit) -fonc(iceigcm_refS)) * dpdiceis(n)
	682	& + (fonc2s(icedsfit) -fonc2s(icedgcm_ref)) * dpdes(n)
	683	pcaps_new= max(pcaps_new,0.)
[2567]	684
[2338]	685	call sol2ls(sol(n),solconv)
	686
	687	write(41,) sol(n), solconv, alphavl1(n)
	688	& ( pfit - pcapn_new - pcaps_new)
	689	write(42,*) sol(n),solconv, pfit - pcapn_new - pcaps_new,
	690	& pcapn_new, pcaps_new, pfit
	691	end do
	692	close(41)
	693	close(42)
	694
	695	end
	696
	697	c *****************************************************************
	698	SUBROUTINE RUNAVE (x,y,nmax,xperiod,xave,xsmooth,ysmooth,nsm)
	699
	700
	701	IMPLICIT NONE
	702
	703	c -----------------------------------------------------------
	704	c Computing runnig average ysmooth on xsmooth coordinate for a periodic
	705	c field y in coordinate x PERIODIC between 0 and xperiod
	706	c Averaging box size : xave
	707	c F.Forget 1999
	708	c -----------------------------------------------------------
	709
	710	integer nmax,nsm
	711	real y(nmax), ysmooth(nsm)
	712	real x(nmax), xsmooth(nsm)
	713	real xperiod,xave
	714	integer n,i, nave, imin
	715
	716	integer nbig
	717	parameter (nbig=99999999)
	718	real xx(nbig)
	719
	720	c -----------------------------------------------------------
	721
	722
[2567]	723	if (nbig.lt.3*nmax) then
	724	write(,) 'Must increase nbig in runave'
	725	stop
	726	endif
[2338]	727
	728	c Reindexation des donnees
	729	do n=1,nmax
	730	xx(n) = x(n) - xperiod
	731	xx(n+nmax) = x(n)
	732	xx(n+2*nmax) = x(n) +xperiod
	733	end do
	734
	735	c Moyenne glissante
	736	imin=1
	737	do n=1,nsm
	738	ysmooth(n) =0.
	739	nave =0
	740	do i=imin,3*nmax
	741	if (xx(i).ge.(xsmooth(n)-0.5*xave)) then
	742	if (xx(i).gt.(xsmooth(n)+0.5*xave)) goto 999
	743	ysmooth(n) = ysmooth(n) + y(mod(i-1,nmax)+1)
	744	nave = nave +1
	745	end if
	746	end do
	747	999 continue
	748	imin = i -1 -nave
	749	ysmooth(n) = ysmooth(n)/ float (nave)
	750	end do
	751	end
	752
	753	subroutine ls2sol(ls,sol)
	754
	755	implicit none
	756	!================================================================
	757	! Arguments:
	758	!================================================================
	759	real,intent(in) :: ls
	760	real,intent(out) :: sol
	761
	762	!================================================================
	763	! Local:
	764	!================================================================
	765	double precision xref,zx0,zteta,zz
	766	!xref: mean anomaly, zteta: true anomaly, zx0: eccentric anomaly
	767	double precision year_day
	768	double precision peri_day,timeperi,e_elips
	769	double precision pi,degrad
	770	parameter (year_day=668.6d0) ! number of sols in a martian year
	771	parameter (peri_day=485.35d0) ! date (in sols) of perihelion
	772	!timeperi: 2pi( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99
	773	parameter (timeperi=1.90258341759902d0)
	774	parameter (e_elips=0.0934d0) ! eccentricity of orbit
	775	parameter (pi=3.14159265358979d0)
	776	parameter (degrad=57.2957795130823d0)
	777
	778	if (abs(ls).lt.1.0e-5) then
	779	if (ls.ge.0.0) then
	780	sol = 0.0
	781	else
	782	sol = real(year_day)
	783	end if
	784	return
	785	end if
	786
	787	zteta = ls/degrad + timeperi
	788	zx0 = 2.0datan(dtan(0.5zteta)/dsqrt((1.+e_elips)/(1.-e_elips)))
	789	xref = zx0-e_elips*dsin(zx0)
	790	zz = xref/(2.*pi)
	791	sol = real(zz*year_day + peri_day)
	792	if (sol.lt.0.0) sol = sol + real(year_day)
	793	if (sol.ge.year_day) sol = sol - real(year_day)
	794
	795
	796	end subroutine ls2sol
	797
	798	subroutine sol2ls(sol,Ls)
	799	!==============================================================================
	800	! Purpose:
	801	! Convert a date/time, given in sol (martian day),
	802	! into solar longitude date/time, in Ls (in degrees),
	803	! where sol=0 is (by definition) the northern hemisphere
	804	! spring equinox (where Ls=0).
	805	!==============================================================================
	806	! Notes:
	807	! Even though "Ls" is cyclic, if "sol" is greater than N (martian) year,
	808	! "Ls" will be increased by N*360
	809	! Won't work as expected if sol is negative (then again,
	810	! why would that ever happen?)
	811	!==============================================================================
	812
	813	implicit none
	814
	815	!==============================================================================
	816	! Arguments:
	817	!==============================================================================
	818	real,intent(in) :: sol
	819	real,intent(out) :: Ls
	820
	821	!==============================================================================
	822	! Local variables:
	823	!==============================================================================
	824	real year_day,peri_day,timeperi,e_elips,twopi,degrad
	825	data year_day /669./ ! # of sols in a martian year
	826	data peri_day /485.0/
	827	data timeperi /1.9082314/
	828	data e_elips /0.093358/
	829	data twopi /6.2831853/ ! 2.*pi
	830	data degrad /57.2957795/ ! pi/180
	831
	832	real zanom,xref,zx0,zdx,zteta,zz
	833
	834	integer count_years
	835	integer iter
	836
	837	!==============================================================================
	838	! 1. Compute Ls
	839	!==============================================================================
	840
	841	zz=(sol-peri_day)/year_day
	842	zanom=twopi*(zz-nint(zz))
	843	xref=abs(zanom)
	844
	845	! The equation zx0 - e * sin (zx0) = xref, solved by Newton
	846	zx0=xref+e_elips*sin(xref)
	847	do iter=1,20 ! typically, 2 or 3 iterations are enough
	848	zdx=-(zx0-e_elipssin(zx0)-xref)/(1.-e_elipscos(zx0))
	849	zx0=zx0+zdx
	850	if(abs(zdx).le.(1.e-7)) then
	851	! write(,)'iter:',iter,' \|zdx\|:',abs(zdx)
	852	exit
	853	endif
	854	enddo
	855
	856	if(zanom.lt.0.) zx0=-zx0
	857
	858	zteta=2.atan(sqrt((1.+e_elips)/(1.-e_elips))tan(zx0/2.))
	859	Ls=zteta-timeperi
	860
	861	if(Ls.lt.0.) then
	862	Ls=Ls+twopi
	863	else
	864	if(Ls.gt.twopi) then
	865	Ls=Ls-twopi
	866	endif
	867	endif
	868
	869	Ls=degrad*Ls
	870	! Ls is now in degrees
	871
	872	!==============================================================================
	873	! 1. Account for (eventual) years included in input date/time sol
	874	!==============================================================================
	875
	876	count_years=0 ! initialize
	877	zz=sol ! use "zz" to store (and work on) the value of sol
	878	do while (zz.ge.year_day)
	879	count_years=count_years+1
	880	zz=zz-year_day
	881	enddo
	882
	883	! Add 360 degrees to Ls for every year
	884	if (count_years.ne.0) then
	885	Ls=Ls+360.*count_years
	886	endif
	887
	888
	889	end subroutine sol2ls
	890

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

Download in other formats:

Original Format