Context Navigation

inifis_mod.F90 @ 3557

Last change on this file since 3557 was 3557, checked in by debatzbr, 2 weeks ago
Miscellaneous cleans + Set-up the physics for the implementation of the microphysical model.
File size: 59.7 KB

Rev	Line
[3184]	1	MODULE inifis_mod
	2	IMPLICIT NONE
	3
	4	CONTAINS
	5
	6	SUBROUTINE inifis(ngrid,nlayer,nq, &
	7	day_ini,pdaysec,nday,ptimestep, &
	8	plat,plon,parea, &
	9	prad,pg,pr,pcpp)
	10
	11	use init_print_control_mod, only: init_print_control
	12	use radinc_h, only: ini_radinc_h, naerkind
	13	use radcommon_h, only: ini_radcommon_h
	14	use radii_mod, only: radfixed, Nmix_n2
[3557]	15	use datafile_mod, only: datadir,config_mufi,hazeprop_file,hazerad_file,hazemmr_file,hazedens_file
[3184]	16	use comdiurn_h, only: sinlat, coslat, sinlon, coslon
	17	use comgeomfi_h, only: totarea, totarea_planet
	18	use comsoil_h, only: ini_comsoil_h, nsoilmx, lay1_soil, alpha_soil
	19	use time_phylmdz_mod, only: ecritphy,day_step,iphysiq, &
	20	init_time, daysec, dtphys
	21	use comcstfi_mod, only: rad, cpp, g, r, rcp, &
	22	mugaz, pi, avocado
	23	use planetwide_mod, only: planetwide_sumval
	24	use callkeys_mod
[3195]	25	use surfdat_h
[3184]	26	use wstats_mod, only: callstats
	27	use ioipsl_getin_p_mod, only : getin_p
	28	use mod_phys_lmdz_para, only : is_parallel, is_master, bcast
	29
	30	!=======================================================================
	31	!
	32	! purpose:
	33	! -------
	34	!
[3197]	35	! Initialisation for the physical parametrisations of the LMD
[3184]	36	! Generic Model.
	37	!
	38	! author: Frederic Hourdin 15 / 10 /93
	39	! -------
	40	! modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
	41	! Ehouarn Millour (oct. 2008) tracers are now identified
	42	! by their names and may not be contiguously
	43	! stored in the q(:,:,:,:) array
	44	! E.M. (june 2009) use getin routine to load parameters
	45	!
	46	!
	47	! arguments:
	48	! ----------
	49	!
	50	! input:
	51	! ------
	52	!
	53	! ngrid Size of the horizontal grid.
	54	! All internal loops are performed on that grid.
	55	! nlayer Number of vertical layers.
	56	! pdayref Day of reference for the simulation
	57	! pday Number of days counted from the North. Spring
	58	! equinoxe.
	59	!
	60	!=======================================================================
	61	!
	62	!-----------------------------------------------------------------------
	63	! declarations:
	64	! -------------
	65	use datafile_mod, only: datadir
	66	use ioipsl_getin_p_mod, only: getin_p
	67	IMPLICIT NONE
	68
	69
	70
	71	REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep
	72	INTEGER,INTENT(IN) :: nday
	73	INTEGER,INTENT(IN) :: ngrid,nlayer,nq
	74	REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid)
	75	integer,intent(in) :: day_ini
	76
	77	INTEGER ig
	78	CHARACTER(len=20) :: rname="inifis" ! routine name, for messages
[3197]	79
[3184]	80	EXTERNAL iniorbit,orbite
	81	EXTERNAL SSUM
	82	REAL SSUM
[3197]	83
[3184]	84	! Initialize flags lunout, prt_level, debug (in print_control_mod)
	85	CALL init_print_control
	86
	87	! initialize constants in comcstfi_mod
	88	rad=prad
	89	cpp=pcpp
	90	g=pg
	91	r=pr
	92	rcp=r/cpp
	93	mugaz=8.314*1000./pr ! dummy init
	94	pi=2.*asin(1.)
	95	avocado = 6.02214179e23 ! added by RW
	96
	97	! Initialize some "temporal and calendar" related variables
	98	#ifndef MESOSCALE
	99	CALL init_time(day_ini,pdaysec,nday,ptimestep)
	100	#endif
	101
	102	! read in some parameters from "run.def" for physics,
	103	! or shared between dynamics and physics.
	104	call getin_p("ecritphy",ecritphy) ! frequency of outputs in physics,
	105	! in dynamical steps
	106	call getin_p("day_step",day_step) ! number of dynamical steps per day
	107	call getin_p("iphysiq",iphysiq) ! call physics every iphysiq dyn step
	108
	109	! do we read a startphy.nc file? (default: .true.)
	110	call getin_p("startphy_file",startphy_file)
[3197]	111
[3184]	112	! --------------------------------------------------------------
	113	! Reading the "callphys.def" file controlling some key options
	114	! --------------------------------------------------------------
[3197]	115
[3184]	116	IF (is_master) THEN
	117	! check if 'callphys.def' file is around
	118	inquire(file="callphys.def",exist=iscallphys)
	119	write(,) trim(rname)//": iscallphys=",iscallphys
	120	ENDIF
	121	call bcast(iscallphys)
	122
	123	IF(iscallphys) THEN
	124
	125	if (is_master) then
	126	write(,)
	127	write(,)'--------------------------------------------'
	128	write(,)trim(rname)//': Parameters for the physics (callphys.def)'
	129	write(,)'--------------------------------------------'
	130	endif
	131
	132	if (is_master) write(,) trim(rname)//&
	133	": Directory where external input files are:"
	134	! default 'datadir' is set in "datadir_mod"
	135	call getin_p("datadir",datadir) ! default path
	136	if (is_master) write(,) trim(rname)//": datadir = ",trim(datadir)
	137
[3193]	138	if (is_master) write(,) "Haze optical properties datafile"
[3329]	139	hazeprop_file="optprop_tholins_fractal050nm.dat" ! default file
	140	call getin_p("hazeprop_file",hazeprop_file)
	141	if (is_master) write(,) trim(rname)//" hazeprop_file = ",trim(hazeprop_file)
[3193]	142
[3329]	143	if (is_master) write(,) "Haze radii datafile"
	144	hazerad_file="hazerad.txt" ! default file
	145	call getin_p("hazerad_file",hazerad_file)
	146	if (is_master) write(,) trim(rname)//" hazerad_file = ",trim(hazerad_file)
[3193]	147
[3329]	148	if (is_master) write(,) "Haze mmr datafile"
[3353]	149	hazemmr_file="None" ! default file
[3329]	150	call getin_p("hazemmr_file",hazemmr_file)
	151	if (is_master) write(,) trim(rname)//" hazemmr_file = ",trim(hazemmr_file)
[3353]	152	if (is_master) write(,) "Haze density datafile"
	153	hazedens_file="None" ! default file
	154	call getin_p("hazedens_file",hazedens_file)
	155	if (is_master) write(,) trim(rname)//" hazedens_file = ",trim(hazedens_file)
[3329]	156
[3184]	157	if (is_master) write(,) trim(rname)//&
	158	": Run with or without tracer transport ?"
	159	tracer=.false. ! default value
	160	call getin_p("tracer",tracer)
	161	if (is_master) write(,) trim(rname)//": tracer = ",tracer
	162
	163	if (is_master) write(,) trim(rname)//&
	164	": Run with or without atm mass update "//&
	165	" due to tracer evaporation/condensation?"
	166	mass_redistrib=.false. ! default value
	167	call getin_p("mass_redistrib",mass_redistrib)
	168	if (is_master) write(,) trim(rname)//&
	169	": mass_redistrib = ",mass_redistrib
	170
	171	if (is_master) then
	172	write(,) trim(rname)//": Diurnal cycle ?"
	173	write(,) trim(rname)//&
	174	": (if diurnal=false, diurnal averaged solar heating)"
	175	endif
	176	diurnal=.true. ! default value
	177	call getin_p("diurnal",diurnal)
	178	if (is_master) write(,) trim(rname)//": diurnal = ",diurnal
	179
	180	if (is_master) then
	181	write(,) trim(rname)//": Seasonal cycle ?"
	182	write(,) trim(rname)//&
	183	": (if season=false, Ls stays constant, to value ", &
	184	"set in 'start'"
	185	endif
	186	season=.true. ! default value
	187	call getin_p("season",season)
	188	if (is_master) write(,) trim(rname)//": season = ",season
[3197]	189
[3228]	190	if (is_master) then
	191	write(,) trim(rname)//": Fast mode (nogcm) ?"
	192	endif
	193	fast=.false. ! default value
	194	call getin_p("fast",fast)
	195	if (is_master) write(,) trim(rname)//": fast = ",fast
	196
[3237]	197	if (is_master) write(,) trim(rname)//"Run with Triton orbit ?"
	198	triton=.false. ! default value
	199	call getin_p("triton",triton)
	200	if (is_master) write(,) trim(rname)//" triton = ",triton
	201	convergeps=.false. ! default value
	202	call getin_p("convergeps",convergeps)
	203	if (is_master) write(,) trim(rname)//" convergeps = ",convergeps
[3539]	204
[3412]	205	conservn2=.false. ! default value
	206	call getin_p("conservn2",conservn2)
	207	if (is_master) write(,) trim(rname)//" conservn2 = ",conservn2
[3539]	208	conservch4=.false. ! default value
	209	call getin_p("conservch4",conservch4)
	210	if (is_master) write(,) trim(rname)//" conservch4 = ",conservch4
[3258]	211
[3237]	212	if (is_master) write(,) trim(rname)//"KBO runs (eris, makemake) ?"
	213	kbo=.false. ! default value
	214	call getin_p("kbo",kbo)
	215	if (is_master) write(,) trim(rname)//" kbo = ",kbo
[3258]	216
[3237]	217	if (is_master) write(,) trim(rname)//"Specific paleo run ?"
	218	paleo=.false. ! default value
	219	call getin_p("paleo",paleo)
	220	if (is_master) write(,) trim(rname)//" paleo = ",paleo
[3258]	221
[3237]	222	if (is_master) write(,) trim(rname)//"paleoclimate step (Earth years) "
	223	paleoyears=10000. ! default value
	224	call getin_p("paleoyears",paleoyears)
	225	if (is_master) write(,) trim(rname)//" paleoyears = ",paleoyears
[3258]	226
[3184]	227	if (is_master) write(,) trim(rname)//&
	228	": No seasonal cycle: initial day to lock the run during restart"
	229	noseason_day=0.0 ! default value
	230	call getin_p("noseason_day",noseason_day)
	231	if (is_master) write(,) trim(rname)//": noseason_day=", noseason_day
	232
	233	if (is_master) write(,) trim(rname)//&
	234	": Write some extra output to the screen ?"
	235	lwrite=.false. ! default value
	236	call getin_p("lwrite",lwrite)
	237	if (is_master) write(,) trim(rname)//": lwrite = ",lwrite
	238
	239	if (is_master) write(,) trim(rname)//&
	240	": Save statistics in file stats.nc ?"
	241	callstats=.true. ! default value
	242	call getin_p("callstats",callstats)
	243	if (is_master) write(,) trim(rname)//": callstats = ",callstats
	244
	245	if (is_master) write(,) trim(rname)//&
	246	": Test energy conservation of model physics ?"
	247	enertest=.false. ! default value
	248	call getin_p("enertest",enertest)
	249	if (is_master) write(,) trim(rname)//": enertest = ",enertest
	250
	251	if (is_master) write(,) trim(rname)//": call radiative transfer ?"
	252	callrad=.true. ! default value
	253	call getin_p("callrad",callrad)
	254	if (is_master) write(,) trim(rname)//": callrad = ",callrad
	255
	256	if (is_master) write(,) trim(rname)//&
	257	": call correlated-k radiative transfer ?"
	258	corrk=.true. ! default value
	259	call getin_p("corrk",corrk)
	260	if (is_master) write(,) trim(rname)//": corrk = ",corrk
	261
	262	if (is_master) write(,) trim(rname)//&
	263	": prohibit calculations outside corrk T grid?"
	264	strictboundcorrk=.true. ! default value
	265	call getin_p("strictboundcorrk",strictboundcorrk)
	266	if (is_master) write(,) trim(rname)//&
	267	": strictboundcorrk = ",strictboundcorrk
	268
	269	if (is_master) write(,) trim(rname)//&
	270	": prohibit calculations outside CIA T grid?"
	271	strictboundcia=.true. ! default value
	272	call getin_p("strictboundcia",strictboundcia)
	273	if (is_master) write(,) trim(rname)//&
	274	": strictboundcia = ",strictboundcia
	275
	276	if (is_master) write(,) trim(rname)//&
	277	": Minimum atmospheric temperature for Planck function integration ?"
	278	tplanckmin=30.0 ! default value
	279	call getin_p("tplanckmin",tplanckmin)
	280	if (is_master) write(,) trim(rname)//": tplanckmin = ",tplanckmin
[3197]	281
[3184]	282	if (is_master) write(,) trim(rname)//&
	283	": Maximum atmospheric temperature for Planck function integration ?"
	284	tplanckmax=1500.0 ! default value
	285	call getin_p("tplanckmax",tplanckmax)
	286	if (is_master) write(,) trim(rname)//": tplanckmax = ",tplanckmax
[3197]	287
[3184]	288	if (is_master) write(,) trim(rname)//&
	289	": Temperature step for Planck function integration ?"
	290	dtplanck=0.1 ! default value
	291	call getin_p("dtplanck",dtplanck)
	292	if (is_master) write(,) trim(rname)//": dtplanck = ",dtplanck
[3197]	293
[3184]	294	if (is_master) write(,) trim(rname)//&
	295	": call gaseous absorption in the visible bands?"//&
	296	" (matters only if callrad=T)"
	297	callgasvis=.false. ! default value
	298	call getin_p("callgasvis",callgasvis)
	299	if (is_master) write(,) trim(rname)//": callgasvis = ",callgasvis
[3197]	300
[3184]	301	if (is_master) write(,) trim(rname)//&
	302	": call continuum opacities in radiative transfer ?"//&
	303	" (matters only if callrad=T)"
	304	continuum=.true. ! default value
	305	call getin_p("continuum",continuum)
	306	if (is_master) write(,) trim(rname)//": continuum = ",continuum
	307
	308	if (is_master) write(,) trim(rname)//&
	309	": call turbulent vertical diffusion ?"
[3539]	310	calldifv=.false. ! default value
[3184]	311	call getin_p("calldifv",calldifv)
	312	if (is_master) write(,) trim(rname)//": calldifv = ",calldifv
	313
	314	if (is_master) write(,) trim(rname)//": use turbdiff instead of vdifc ?"
[3539]	315	UseTurbDiff=.false. ! default value
[3184]	316	call getin_p("UseTurbDiff",UseTurbDiff)
	317	if (is_master) write(,) trim(rname)//": UseTurbDiff = ",UseTurbDiff
	318
	319	if (is_master) write(,) trim(rname)//": call convective adjustment ?"
[3539]	320	calladj=.false. ! default value
[3184]	321	call getin_p("calladj",calladj)
	322	if (is_master) write(,) trim(rname)//": calladj = ",calladj
	323
	324	if (is_master) write(,) trim(rname)//&
	325	": Radiative timescale for Newtonian cooling (in Earth days)?"
	326	tau_relax=30. ! default value
	327	call getin_p("tau_relax",tau_relax)
	328	if (is_master) write(,) trim(rname)//": tau_relax = ",tau_relax
	329	tau_relax=tau_relax243600 ! convert Earth days --> seconds
	330
	331	if (is_master) write(,)trim(rname)//&
	332	": call thermal conduction in the soil ?"
	333	callsoil=.true. ! default value
	334	call getin_p("callsoil",callsoil)
	335	if (is_master) write(,) trim(rname)//": callsoil = ",callsoil
[3197]	336
[3184]	337	if (is_master) write(,)trim(rname)//&
	338	": Rad transfer is computed every iradia", &
	339	" physical timestep"
	340	iradia=1 ! default value
	341	call getin_p("iradia",iradia)
	342	if (is_master) write(,)trim(rname)//": iradia = ",iradia
[3197]	343
[3184]	344	if (is_master) write(,)trim(rname)//": Rayleigh scattering ?"
	345	rayleigh=.false.
	346	call getin_p("rayleigh",rayleigh)
	347	if (is_master) write(,)trim(rname)//": rayleigh = ",rayleigh
	348
	349	if (is_master) write(,) trim(rname)//&
	350	": Use blackbody for stellar spectrum ?"
	351	stelbbody=.false. ! default value
	352	call getin_p("stelbbody",stelbbody)
	353	if (is_master) write(,) trim(rname)//": stelbbody = ",stelbbody
	354
	355	if (is_master) write(,) trim(rname)//": Stellar blackbody temperature ?"
	356	stelTbb=5800.0 ! default value
	357	call getin_p("stelTbb",stelTbb)
	358	if (is_master) write(,) trim(rname)//": stelTbb = ",stelTbb
	359
	360	if (is_master) write(,) trim(rname)//": Output mean OLR in 1D?"
	361	meanOLR=.false.
	362	call getin_p("meanOLR",meanOLR)
	363	if (is_master) write(,) trim(rname)//": meanOLR = ",meanOLR
	364
	365	if (is_master) write(,)trim(rname)//": Output spectral OLR in 3D?"
	366	specOLR=.false.
	367	call getin_p("specOLR",specOLR)
	368	if (is_master) write(,)trim(rname)//": specOLR = ",specOLR
	369
	370	if (is_master) write(,)trim(rname)//&
	371	": Output diagnostic optical thickness attenuation exp(-tau)?"
	372	diagdtau=.false.
	373	call getin_p("diagdtau",diagdtau)
	374	if (is_master) write(,)trim(rname)//": diagdtau = ",diagdtau
	375
	376	if (is_master) write(,)trim(rname)//": Operate in kastprof mode?"
	377	kastprof=.false.
	378	call getin_p("kastprof",kastprof)
	379	if (is_master) write(,)trim(rname)//": kastprof = ",kastprof
	380
	381	if (is_master) write(,)trim(rname)//&
	382	": Uniform absorption in radiative transfer?"
	383	graybody=.false.
	384	call getin_p("graybody",graybody)
	385	if (is_master) write(,)trim(rname)//": graybody = ",graybody
	386
	387	! Soil model
	388	if (is_master) write(,)trim(rname)//&
	389	": Number of sub-surface layers for soil scheme?"
	390	! default value of nsoilmx set in comsoil_h
	391	call getin_p("nsoilmx",nsoilmx)
	392	if (is_master) write(,)trim(rname)//": nsoilmx=",nsoilmx
[3197]	393
[3184]	394	if (is_master) write(,)trim(rname)//&
	395	": Thickness of topmost soil layer (m)?"
	396	! default value of lay1_soil set in comsoil_h
	397	call getin_p("lay1_soil",lay1_soil)
	398	if (is_master) write(,)trim(rname)//": lay1_soil = ",lay1_soil
[3197]	399
[3184]	400	if (is_master) write(,)trim(rname)//&
	401	": Coefficient for soil layer thickness distribution?"
	402	! default value of alpha_soil set in comsoil_h
	403	call getin_p("alpha_soil",alpha_soil)
	404	if (is_master) write(,)trim(rname)//": alpha_soil = ",alpha_soil
	405
	406	! Chemistry in the thermosphere
	407	if (is_master) write(,) trim(rname)//": Use deposition ?"
	408	depos=.false. ! default value
	409	call getin_p("depos",depos)
	410	if (is_master) write(,) trim(rname)//": depos = ",depos
	411
	412	if (is_master) write(,)trim(rname)//": Production of haze ?"
	413	haze=.false. ! default value
	414	call getin_p("haze",haze)
	415	if (is_master) write(,)trim(rname)//": haze = ",haze
	416
[3455]	417
	418
	419	if (is_master) write(,)trim(rname)// "call thermal conduction ?"
	420	callconduct=.false. ! default value
	421	call getin_p("callconduct",callconduct)
	422	if (is_master) write(,)trim(rname)// " callconduct = ",callconduct
	423
	424	if (is_master) write(,)trim(rname)// "call phitop ?"
	425	phitop=0. ! default value
	426	call getin_p("phitop",phitop)
	427	if (is_master) write(,)trim(rname)// " phitop = ",phitop
	428
	429	if (is_master) write(,)trim(rname)// "call molecular viscosity ?"
	430	callmolvis=.false. ! default value
	431	call getin_p("callmolvis",callmolvis)
	432	if (is_master) write(,)trim(rname)// " callmolvis = ",callmolvis
	433
	434	if (is_master) write(,)trim(rname)// "call molecular diffusion ?"
	435	callmoldiff=.false. ! default value
	436	call getin_p("callmoldiff",callmoldiff)
	437	if (is_master) write(,)trim(rname)// " callmoldiff = ",callmoldiff
	438
[3184]	439	! Global1D mean and solar zenith angle
	440
	441	if (ngrid.eq.1) then
	442	PRINT*, 'Simulate global averaged conditions ?'
	443	global1d = .false. ! default value
	444	call getin_p("global1d",global1d)
	445	write(,) "global1d = ",global1d
[3197]	446
[3184]	447	! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
	448	if (global1d.and.diurnal) then
	449	call abort_physic(rname,'if global1d is true, diurnal must be set to false',1)
	450	endif
	451
	452	if (global1d) then
	453	PRINT *,'Solar Zenith angle (deg.) ?'
	454	PRINT *,'(assumed for averaged solar flux S/4)'
	455	szangle=60.0 ! default value
	456	call getin_p("szangle",szangle)
	457	write(,) "szangle = ",szangle
	458	endif
	459	endif ! of if (ngrid.eq.1)
	460
	461	! Test of incompatibility:
	462	! if kastprof used, we must be in 1D
	463	if (kastprof.and.(ngrid.gt.1)) then
	464	call abort_physic(rname,'kastprof can only be used in 1D!',1)
	465	endif
	466
	467	if (is_master) write(,)trim(rname)//&
	468	": Stratospheric temperature for kastprof mode?"
	469	Tstrat=167.0
	470	call getin_p("Tstrat",Tstrat)
	471	if (is_master) write(,)trim(rname)//": Tstrat = ",Tstrat
	472
	473	if (is_master) write(,)trim(rname)//": Remove lower boundary?"
	474	nosurf=.false.
	475	call getin_p("nosurf",nosurf)
	476	if (is_master) write(,)trim(rname)//": nosurf = ",nosurf
	477
	478	! Tests of incompatibility:
	479	if (nosurf.and.callsoil) then
	480	if (is_master) then
	481	write(,)trim(rname)//'nosurf not compatible with soil scheme!'
	482	write(,)trim(rname)//'... got to make a choice!'
	483	endif
	484	call abort_physic(rname,"nosurf not compatible with soil scheme!",1)
	485	endif
	486
	487	if (is_master) write(,)trim(rname)//": Add an internal heat flux?", &
	488	"... matters only if callsoil=F"
	489	intheat=0.
	490	call getin_p("intheat",intheat)
	491	if (is_master) write(,)trim(rname)//": intheat = ",intheat
	492
	493
	494	if (is_master) write(,)trim(rname)//": Test physics timescale in 1D?"
	495	testradtimes=.false.
	496	call getin_p("testradtimes",testradtimes)
	497	if (is_master) write(,)trim(rname)//": testradtimes = ",testradtimes
	498
	499	! Test of incompatibility:
	500	! if testradtimes used, we must be in 1D
	501	if (testradtimes.and.(ngrid.gt.1)) then
	502	call abort_physic(rname,'testradtimes can only be used in 1D!',1)
	503	endif
	504
	505	if (is_master) write(,)trim(rname)//": Default planetary temperature?"
	506	tplanet=215.0
	507	call getin_p("tplanet",tplanet)
	508	if (is_master) write(,)trim(rname)//": tplanet = ",tplanet
	509
	510	if (is_master) write(,)trim(rname)//": Which star?"
	511	startype=1 ! default value = Sol
	512	call getin_p("startype",startype)
	513	if (is_master) write(,)trim(rname)//": startype = ",startype
	514
	515	if (is_master) write(,)trim(rname)//": Value of stellar flux at 1 AU?"
	516	Fat1AU=1356.0 ! default value = Sol today
	517	call getin_p("Fat1AU",Fat1AU)
	518	if (is_master) write(,)trim(rname)//": Fat1AU = ",Fat1AU
	519
	520
	521	! TRACERS:
	522
	523	if (is_master) write(,)trim(rname)//": fixed radii for Cloud particles?"
	524	radfixed=.false. ! default value
	525	call getin_p("radfixed",radfixed)
	526	if (is_master) write(,)trim(rname)//": radfixed = ",radfixed
	527
	528	if(kastprof)then
	529	radfixed=.true.
[3197]	530	endif
[3184]	531
[3195]	532	if (is_master) write(,)trim(rname)//&
[3184]	533	"Number of radiatively active aerosols:"
	534	naerkind=0 ! default value
	535	call getin_p("naerkind",naerkind)
	536	if (is_master) write(,)trim(rname)//": naerkind = ",naerkind
	537
[3197]	538
[3195]	539	!***************************************************************
	540	!! TRACERS options
[3184]	541
[3195]	542	if (is_master)write(,)trim(rname)//&
	543	"call N2 condensation ?"
[3197]	544	n2cond=.true. ! default value
	545	call getin_p("n2cond",n2cond)
[3195]	546	if (is_master)write(,)trim(rname)//&
[3197]	547	" n2cond = ",n2cond
[3184]	548
[3195]	549	if (is_master)write(,)trim(rname)//&
	550	"call N2 cloud condensation ?"
	551	condensn2=.false. ! default value
	552	call getin_p("condensn2",condensn2)
	553	if (is_master)write(,)trim(rname)//&
	554	"condensn2 = ",condensn2
[3184]	555
[3195]	556	if (is_master)write(,)trim(rname)//&
	557	"call no N2 frost formation?"
	558	no_n2frost=.false. ! default value
	559	call getin_p("no_n2frost",no_n2frost)
	560	if (is_master)write(,)trim(rname)//&
	561	"no_n2frost = ",no_n2frost
[3184]	562
[3195]	563	if (is_master)write(,)trim(rname)//&
	564	"N2 condensation subtimestep?"
	565	nbsub=20 ! default value
	566	call getin_p("nbsub",nbsub)
	567	if (is_master)write(,)trim(rname)//&
	568	" nbsub = ",nbsub
[3184]	569
[3195]	570	if (is_master)write(,)trim(rname)//&
	571	"Gravitationnal sedimentation ?"
	572	sedimentation=.true. ! default value
	573	call getin_p("sedimentation",sedimentation)
	574	if (is_master)write(,)trim(rname)//&
	575	" sedimentation = ",sedimentation
[3184]	576
[3195]	577	if (is_master)write(,)trim(rname)//&
	578	"Compute glacial flow ?"
	579	glaflow=.false. ! default value
	580	call getin_p("glaflow",glaflow)
	581	if (is_master)write(,)trim(rname)//&
	582	" glaflow = ",glaflow
	583
	584	if (is_master)write(,)trim(rname)//&
	585	"Compute methane cycle ?"
	586	methane=.false. ! default value
	587	call getin_p("methane",methane)
	588	if (is_master)write(,)trim(rname)//&
	589	" methane = ",methane
	590	condmetsurf=.true. ! default value
	591	call getin_p("condmetsurf",condmetsurf)
	592	if (is_master)write(,)trim(rname)//&
	593	" condmetsurf = ",condmetsurf
	594
	595	if (is_master)write(,)trim(rname)//&
	596	"Compute methane clouds ?"
	597	metcloud=.false. ! default value
	598	call getin_p("metcloud",metcloud)
	599	if (is_master)write(,)trim(rname)//&
	600	" metcloud = ",metcloud
	601
	602	if (is_master)write(,)trim(rname)//&
	603	"Compute CO cycle ?"
	604	carbox=.false. ! default value
	605	call getin_p("carbox",carbox)
	606	if (is_master)write(,)trim(rname)//&
	607	" carbox = ",carbox
	608	condcosurf=.true. ! default value
	609	call getin_p("condcosurf",condcosurf)
	610	if (is_master)write(,)trim(rname)//&
	611	" condcosurf = ",condcosurf
	612
	613	if (is_master)write(,)trim(rname)//&
	614	"Compute CO clouds ?"
	615	monoxcloud=.false. ! default value
	616	call getin_p("monoxcloud",monoxcloud)
	617	if (is_master)write(,)trim(rname)//&
	618	" monoxcloud = ",monoxcloud
	619
	620	if (is_master)write(,)trim(rname)//&
	621	"atmospheric redistribution (s):"
	622	tau_n2=1. ! default value
	623	call getin_p("tau_n2",tau_n2)
	624	if (is_master)write(,)trim(rname)//&
	625	" tau_n2 = ",tau_n2
	626	tau_ch4=1.E7 ! default value
	627	call getin_p("tau_ch4",tau_ch4)
	628	if (is_master)write(,)trim(rname)//&
	629	" tau_ch4 = ",tau_ch4
	630	tau_co=1. ! default value
	631	call getin_p("tau_co",tau_co)
	632	if (is_master)write(,)trim(rname)//&
	633	" tau_co = ",tau_co
	634	tau_prechaze=1. ! default value
	635	call getin_p("tau_prechaze",tau_prechaze)
	636	if (is_master)write(,)trim(rname)//&
	637	" tau_prechaze = ",tau_prechaze
	638
	639	if (is_master)write(,)trim(rname)//&
	640	"Day fraction for limited cold trap in SP?"
	641	dayfrac=0. ! default value
	642	call getin_p("dayfrac",dayfrac)
	643	if (is_master)write(,)trim(rname)//&
	644	" dayfrac = ",dayfrac
	645	thresh_non2=0. ! default value
	646	call getin_p("thresh_non2",thresh_non2)
	647	if (is_master)write(,)trim(rname)//&
	648	" thresh_non2 = ",thresh_non2
	649
[3353]	650	! use Pluto.old routines
	651	if (is_master) write(,) trim(rname)//": use vdifc from old Pluto ?"
[3378]	652	oldplutovdifc=.false. ! default value
[3353]	653	call getin_p("oldplutovdifc",oldplutovdifc)
	654	if (is_master) write(,) trim(rname)//": oldplutovdifc = ",oldplutovdifc
	655
	656	if (is_master) write(,) trim(rname)//&
	657	": call pluto.old correlated-k radiative transfer ?"
[3378]	658	oldplutocorrk=.false. ! default value
[3353]	659	call getin_p("oldplutocorrk",oldplutocorrk)
	660	if (is_master) write(,) trim(rname)//": oldplutocorrk = ",oldplutocorrk
	661
	662	if (is_master) write(,) trim(rname)//&
	663	": call pluto.old sedimentation ?"
[3378]	664	oldplutosedim=.false. ! default value
[3353]	665	call getin_p("oldplutosedim",oldplutosedim)
	666	if (is_master) write(,) trim(rname)//": oldplutosedim = ",oldplutosedim
	667
[3195]	668	!***************************************************************
	669	!! Haze options
	670
[3557]	671	! Microphysical moment model
	672	! ~~~~~~~~~~~~~~~~~~~~~~~~~~
	673	if (is_master) write(,) "Run with or without microphysics?"
	674	callmufi=.false. ! default value
	675	call getin_p("callmufi",callmufi)
	676	if (is_master) write(,)" callmufi = ",callmufi
	677
	678	! sanity check
	679	if (callmufi.and.(.not.tracer)) then
	680	print*,"You are running microphysics without tracer"
	681	print*,"Please start again with tracer =.true."
	682	stop
	683	endif
	684
	685	if (is_master) write(,) "Path to microphysical config file?"
	686	config_mufi='datagcm/microphysics/config.cfg' ! default value
	687	call getin_p("config_mufi",config_mufi)
	688	if (is_master) write(,)" config_mufi = ",config_mufi
	689
	690	if (is_master) write(,) "Use haze production from CH4 photolysis or production rate?"
	691	call_haze_prod_pCH4=.false. ! default value
	692	call getin_p("call_haze_prod_pCH4",call_haze_prod_pCH4)
	693	if (is_master) write(,)" call_haze_prod_pCH4 = ",call_haze_prod_pCH4
	694
	695	if (is_master) write(,) "Pressure level of aerosols production (Pa)?"
	696	haze_p_prod=1.0e-2 ! default value
	697	call getin_p("haze_p_prod",haze_p_prod)
	698	if (is_master) write(,)" haze_p_prod = ",haze_p_prod
	699
	700	if (is_master) write(,) "Aerosol production rate (kg.m-2.s-1)?"
	701	haze_tx_prod=9.8e-14 ! default value
	702	call getin_p("haze_tx_prod",haze_tx_prod)
	703	if (is_master) write(,)" haze_tx_prod = ",haze_tx_prod
	704
	705	if (is_master) write(,) "Equivalent radius production (m)?"
	706	haze_rc_prod=1.0e-9 ! default value
	707	call getin_p("haze_rc_prod",haze_rc_prod)
	708	if (is_master) write(,)" haze_rc_prod = ",haze_rc_prod
	709
	710	if (is_master) write(,) "Monomer radius (m)?"
	711	haze_rm=1.0e-8 ! default value
	712	call getin_p("haze_rm",haze_rm)
	713	if (is_master) write(,)" haze_rm = ",haze_rm
	714
	715	if (is_master) write(,) "Aerosol's fractal dimension?"
	716	haze_df=2.0 ! default value
	717	call getin_p("haze_df",haze_df)
	718	if (is_master) write(,)" haze_df = ",haze_df
	719
	720	if (is_master) write(,) "Aerosol density (kg.m-3)?"
	721	haze_rho=800.0 ! default value
	722	call getin_p("haze_rho",haze_rho)
	723	if (is_master) write(,)" haze_rho = ",haze_rho
	724
	725	if (is_master) write(,) "Radius of air molecule (m)?"
	726	air_rad=1.75e-10 ! default value
	727	call getin_p("air_rad",air_rad)
	728	if (is_master) write(,)" air_rad = ",air_rad
	729
	730	! Pluto haze model
	731	! ~~~~~~~~~~~~~~~~
[3195]	732	if (is_master)write(,)trim(rname)//&
	733	"Production of haze ?"
	734	haze=.false. ! default value
	735	call getin_p("haze",haze)
	736	if (is_master)write(,)trim(rname)//&
	737	" haze = ",haze
	738	hazeconservch4=.false. ! conservch4, by default value ch4 is photolyzed
	739	call getin_p("hazeconservch4",hazeconservch4)
	740	if (is_master)write(,)trim(rname)//&
	741	"hazconservch4 = ",hazeconservch4
	742	if (is_master)write(,)trim(rname)//&
	743	"Production of haze (fast version) ?"
	744	fasthaze=.false. ! default value
	745	call getin_p("fasthaze",fasthaze)
	746	if (is_master)write(,)trim(rname)//&
	747	"fasthaze = ",fasthaze
	748
	749	if (is_master)write(,)trim(rname)//&
	750	"Add source of haze ?"
	751	source_haze=.false. ! default value
	752	call getin_p("source_haze",source_haze)
	753	if (is_master)write(,)trim(rname)//&
	754	" source_haze = ",source_haze
	755	mode_hs=0 ! mode haze source default value
	756	call getin_p("mode_hs",mode_hs)
	757	if (is_master)write(,)trim(rname)//&
	758	" mode_hs",mode_hs
	759	kfix=1 ! default value
	760	call getin_p("kfix",kfix)
	761	if (is_master)write(,)trim(rname)//&
	762	"levels kfix",kfix
	763	fracsource=0.1 ! default value
	764	call getin_p("fracsource",fracsource)
	765	if (is_master)write(,)trim(rname)//&
	766	" fracsource",fracsource
	767	latsource=30. ! default value
	768	call getin_p("latsource",latsource)
	769	if (is_master)write(,)trim(rname)//&
	770	" latsource",latsource
	771	lonsource=180. ! default value
	772	call getin_p("lonsource",lonsource)
	773	if (is_master)write(,)trim(rname)//&
	774	" lonsource",lonsource
	775
	776	if (is_master)write(,)trim(rname)//&
	777	"Radiatively active haze ?"
	778	aerohaze=.false. ! default value
	779	call getin_p("aerohaze",aerohaze)
	780	if (is_master)write(,)trim(rname)//&
	781	"aerohaze = ",aerohaze
	782
	783	if (is_master)write(,)trim(rname)//&
	784	"Haze monomer radius ?"
	785	rad_haze=10.e-9 ! default value
	786	call getin_p("rad_haze",rad_haze)
	787	if (is_master)write(,)trim(rname)//&
	788	"rad_haze = ",rad_haze
	789
	790	if (is_master)write(,)trim(rname)//&
	791	"fractal particle ?"
	792	fractal=.false. ! default value
	793	call getin_p("fractal",fractal)
	794	if (is_master)write(,)trim(rname)//&
	795	"fractal = ",fractal
	796	nb_monomer=10 ! default value
	797	call getin_p("nb_monomer",nb_monomer)
	798	if (is_master)write(,)trim(rname)//&
	799	"nb_monomer = ",nb_monomer
[3197]	800
[3195]	801	if (is_master)write(,)trim(rname)//&
[3353]	802	"fixed gaseous CH4 mixing ratio for rad transfer ?"
	803	ch4fix=.false. ! default value
	804	call getin_p("ch4fix",ch4fix)
	805	if (is_master)write(,)trim(rname)//&
	806	" ch4fix = ",ch4fix
	807	if (is_master)write(,)trim(rname)//&
	808	"fixed gaseous CH4 mixing ratio for rad transfer ?"
	809	vmrch4fix=0.5 ! default value
	810	call getin_p("vmrch4fix",vmrch4fix)
	811	if (is_master)write(,)trim(rname)//&
	812	" vmrch4fix = ",vmrch4fix
	813	vmrch4_proffix=.false. ! default value
	814	call getin_p("vmrch4_proffix",vmrch4_proffix)
	815	if (is_master)write(,)trim(rname)//&
	816	" vmrch4_proffix = ",vmrch4_proffix
	817
	818	if (is_master)write(,)trim(rname)//&
	819	"call specific cooling for radiative transfer ?"
	820	cooling=.false. ! default value
	821	call getin_p("cooling",cooling)
	822	if (is_master)write(,)trim(rname)//&
	823	" cooling = ",cooling
	824
	825	if (is_master)write(,)trim(rname)//&
	826	"NLTE correction for methane heating rates?"
	827	nlte=.false. ! default value
	828	call getin_p("nlte",nlte)
	829	if (is_master)write(,)trim(rname)//&
	830	" nlte = ",nlte
	831	strobel=.false. ! default value
	832	call getin_p("strobel",strobel)
	833	if (is_master)write(,)trim(rname)//&
	834	" strobel = ",strobel
	835
	836	if (is_master)write(,)trim(rname)//&
[3195]	837	"fixed radius profile from txt file ?"
	838	haze_radproffix=.false. ! default value
	839	call getin_p("haze_radproffix",haze_radproffix)
	840	if (is_master)write(,)trim(rname)//&
	841	"haze_radproffix = ",haze_radproffix
	842	if (is_master)write(,)trim(rname)//&
	843	"fixed MMR profile from txt file ?"
	844	haze_proffix=.false. ! default value
	845	call getin_p("haze_proffix",haze_proffix)
	846	if (is_master)write(,)trim(rname)//&
	847	"haze_proffix = ",haze_proffix
	848
	849	if (is_master)write(,)trim(rname)//&
	850	"Number mix ratio of haze particles for co clouds:"
	851	Nmix_co=100000. ! default value
	852	call getin_p("Nmix_co",Nmix_co)
	853	if (is_master)write(,)trim(rname)//&
	854	" Nmix_co = ",Nmix_co
	855
	856	if (is_master)write(,)trim(rname)//&
	857	"Number mix ratio of haze particles for ch4 clouds:"
	858	Nmix_ch4=100000. ! default value
	859	call getin_p("Nmix_ch4",Nmix_ch4)
	860	if (is_master)write(,)trim(rname)//&
	861	" Nmix_ch4 = ",Nmix_ch4
	862
	863	!***************************************************************
	864	!! Surface properties
	865
	866	!********* N2 *******************************
	867
	868	if (is_master)write(,)trim(rname)//&
	869	"Mode for changing N2 albedo"
	870	mode_n2=0 ! default value
	871	call getin_p("mode_n2",mode_n2)
	872	if (is_master)write(,)trim(rname)//&
	873	" mode_n2 = ",mode_n2
	874	thres_n2ice=1. ! default value
	875	call getin_p("thres_n2ice",thres_n2ice)
	876	if (is_master)write(,)trim(rname)//&
	877	" thres_n2ice = ",thres_n2ice
	878
	879	if (is_master)write(,)trim(rname)//&
	880	"Diff of N2 albedo with thickness"
	881	deltab=0. ! default value
	882	call getin_p("deltab",deltab)
	883	if (is_master)write(,)trim(rname)//&
	884	" deltab = ",deltab
	885
	886	if (is_master)write(,)trim(rname)//&
	887	"albedo N2 beta "
	888	alb_n2b=0.7 ! default value
	889	call getin_p("alb_n2b",alb_n2b)
	890	if (is_master)write(,)trim(rname)//&
	891	" alb_n2b = ",alb_n2b
	892
	893	if (is_master)write(,)trim(rname)//&
	894	"albedo N2 alpha "
	895	alb_n2a=0.7 ! default value
	896	call getin_p("alb_n2a",alb_n2a)
	897	if (is_master)write(,)trim(rname)//&
	898	" alb_n2a = ",alb_n2a
	899
	900	if (is_master)write(,)trim(rname)//&
	901	"emis N2 beta "
	902	emis_n2b=0.7 ! default value
	903	call getin_p("emis_n2b",emis_n2b)
	904	if (is_master)write(,)trim(rname)//&
	905	" emis_n2b = ",emis_n2b
	906
	907	if (is_master)write(,)trim(rname)//&
	908	"emis N2 alpha "
	909	emis_n2a=0.7 ! default value
	910	call getin_p("emis_n2a",emis_n2a)
	911	if (is_master)write(,)trim(rname)//&
	912	" emis_n2a = ",emis_n2a
	913
	914	!********* CH4 *******************************
	915
	916	if (is_master)write(,)trim(rname)//&
	917	"Mode for changing CH4 albedo"
	918	mode_ch4=0 ! default value
	919	call getin_p("mode_ch4",mode_ch4)
	920	if (is_master)write(,)trim(rname)//&
	921	" mode_ch4 = ",mode_ch4
	922	feedback_met=0 ! default value
	923	call getin_p("feedback_met",feedback_met)
	924	if (is_master)write(,)trim(rname)//&
	925	" feedback_met = ",feedback_met
	926	thres_ch4ice=1. ! default value, mm
	927	call getin_p("thres_ch4ice",thres_ch4ice)
	928	if (is_master)write(,)trim(rname)//&
	929	" thres_ch4ice = ",thres_ch4ice
[3197]	930	fdch4_latn=-1.
	931	fdch4_lats=0.
	932	fdch4_lone=0.
	933	fdch4_lonw=-1.
	934	fdch4_depalb=0.5
	935	fdch4_finalb=0.95
	936	fdch4_maxalb=0.99
	937	fdch4_ampl=200.
	938	fdch4_maxice=100.
[3372]	939	call getin_p("fdch4_latn",fdch4_latn)
	940	call getin_p("fdch4_lats",fdch4_lats)
	941	call getin_p("fdch4_lone",fdch4_lone)
	942	call getin_p("fdch4_lonw",fdch4_lonw)
	943	call getin_p("fdch4_depalb",fdch4_depalb)
	944	call getin_p("fdch4_finalb",fdch4_finalb)
	945	call getin_p("fdch4_maxalb",fdch4_maxalb)
	946	call getin_p("fdch4_maxice",fdch4_maxice)
	947	call getin_p("fdch4_ampl",fdch4_ampl)
[3195]	948	if (is_master)write(,)trim(rname)//&
	949	" Values for albedo feedback = ",fdch4_latn,&
	950	fdch4_lats,fdch4_lone,fdch4_lonw,fdch4_depalb,&
	951	fdch4_finalb,fdch4_maxalb,fdch4_maxice,fdch4_ampl
	952
	953	if (is_master)write(,)trim(rname)//&
	954	"Latitude for diff albedo methane"
	955	metlateq=25. ! default value
	956	call getin_p("metlateq",metlateq)
	957	if (is_master)write(,)trim(rname)//&
	958	" metlateq = ",metlateq
	959
	960	if (is_master)write(,)trim(rname)//&
	961	"Ls1 and Ls2 for change of ch4 albedo"
	962	metls1=-1. ! default value
	963	metls2=-2. ! default value
	964	call getin_p("metls1",metls1)
	965	call getin_p("metls2",metls2)
	966
	967	if (is_master)write(,)trim(rname)//&
	968	"albedo CH4 "
	969	alb_ch4=0.5 ! default value
	970	call getin_p("alb_ch4",alb_ch4)
	971	if (is_master)write(,)trim(rname)//&
	972	" alb_ch4 = ",alb_ch4
	973
	974	if (is_master)write(,)trim(rname)//&
	975	"albedo equatorial CH4 "
	976	alb_ch4_eq=alb_ch4 ! default value
	977	call getin_p("alb_ch4_eq",alb_ch4_eq)
	978	if (is_master)write(,)trim(rname)//&
	979	" alb_ch4_eq = ",alb_ch4_eq
	980
	981	if (is_master)write(,)trim(rname)//&
	982	"albedo south hemis CH4 "
	983	alb_ch4_s=alb_ch4 ! default value
	984	call getin_p("alb_ch4_s",alb_ch4_s)
	985	if (is_master)write(,)trim(rname)//&
	986	" alb_ch4_s = ",alb_ch4_s
	987
	988	if (is_master)write(,)trim(rname)//&
	989	"emis CH4 "
	990	emis_ch4=0.5 ! default value
	991	call getin_p("emis_ch4",emis_ch4)
	992	if (is_master)write(,)trim(rname)//&
	993	" emis_ch4 = ",emis_ch4
	994
	995	if (is_master)write(,)trim(rname)//&
	996	"CH4 lag for n2 sublimation limitation"
	997	ch4lag=.false. ! default value
	998	latlag=-90. ! default value
	999	vmrlag=1. ! default value
	1000	call getin_p("ch4lag",ch4lag)
	1001	call getin_p("latlag",latlag)
	1002	call getin_p("vmrlag",vmrlag)
	1003	if (is_master)write(,)trim(rname)//&
	1004	" ch4lag = ",ch4lag
	1005	if (is_master)write(,)trim(rname)//&
	1006	" latlag = ",latlag
	1007	if (is_master)write(,)trim(rname)//&
	1008	" vmrlag = ",vmrlag
	1009
	1010	if (is_master)write(,)trim(rname)//&
	1011	"Max temperature for surface ?"
	1012	tsurfmax=.false. ! default value
	1013	albmin_ch4=0.3 ! default value
	1014	call getin_p("tsurfmax",tsurfmax)
	1015	call getin_p("albmin_ch4",albmin_ch4)
	1016	if (is_master)write(,)trim(rname)//&
	1017	" tsurfmax = ",tsurfmax
	1018	if (is_master)write(,)trim(rname)//&
	1019	" albmin_ch4 = ",albmin_ch4
	1020
[3275]	1021
	1022	if (is_master)write(,)trim(rname)//&
	1023	"fixed gaseous CH4 mixing ratio for rad transfer ?"
	1024	ch4fix=.false. ! default value
	1025	call getin_p("ch4fix",ch4fix)
	1026	if (is_master)write(,)trim(rname)//&
	1027	" ch4fix = ",ch4fix
	1028	if (is_master)write(,)trim(rname)//&
	1029	"fixed gaseous CH4 mixing ratio for rad transfer ?"
	1030	vmrch4fix=0.5 ! default value
	1031	call getin_p("vmrch4fix",vmrch4fix)
	1032	if (is_master)write(,)trim(rname)//&
	1033	" vmrch4fix = ",vmrch4fix
	1034	vmrch4_proffix=.false. ! default value
	1035	call getin_p("vmrch4_proffix",vmrch4_proffix)
	1036	if (is_master)write(,)trim(rname)//&
	1037	" vmrch4_proffix = ",vmrch4_proffix
	1038
	1039
[3195]	1040	!********* CO *******************************
	1041
	1042	if (is_master)write(,)trim(rname)//&
	1043	"albedo CO "
	1044	alb_co=0.4 ! default value
	1045	call getin_p("alb_co",alb_co)
	1046	if (is_master)write(,)trim(rname)//&
	1047	" alb_co = ",alb_co
	1048	thres_coice=1. ! default value, mm
	1049	call getin_p("thres_coice",thres_coice)
	1050	if (is_master)write(,)trim(rname)//&
	1051	" thres_coice = ",thres_coice
	1052
	1053	if (is_master)write(,)trim(rname)//&
	1054	"emis CO "
	1055	emis_co=0.4 ! default value
	1056	call getin_p("emis_co",emis_co)
	1057	if (is_master)write(,)trim(rname)//&
	1058	" emis_co = ",emis_co
	1059
	1060	!********* THOLINS *******************************
	1061	if (is_master)write(,)trim(rname)//&
	1062	"Mode for changing tholins albedo/emis"
	1063	mode_tholins=0 ! default value
	1064	call getin_p("mode_tholins",mode_tholins)
	1065	if (is_master)write(,)trim(rname)//&
	1066	" mode_tholins = ",mode_tholins
	1067
	1068	if (is_master)write(,)trim(rname)//&
	1069	"albedo tho "
	1070	alb_tho=0.1 ! default value
	1071	call getin_p("alb_tho",alb_tho)
	1072	if (is_master)write(,)trim(rname)//&
	1073	" alb_tho = ",alb_tho
	1074
	1075	if (is_master)write(,)trim(rname)//&
	1076	"albedo tho eq"
	1077	alb_tho_eq=0.1 ! default value
	1078	call getin_p("alb_tho_eq",alb_tho_eq)
	1079	if (is_master)write(,)trim(rname)//&
	1080	" alb_tho_eq = ",alb_tho_eq
	1081
	1082	if (is_master)write(,)trim(rname)//&
	1083	"emis tholins "
	1084	emis_tho=1. ! default value
	1085	call getin_p("emis_tho",emis_tho)
	1086	if (is_master)write(,)trim(rname)//&
	1087	" emis_tho = ",emis_tho
	1088
	1089	if (is_master)write(,)trim(rname)//&
	1090	"emis tholins eq"
	1091	emis_tho_eq=1. ! default value
	1092	call getin_p("emis_tho_eq",emis_tho_eq)
	1093	if (is_master)write(,)trim(rname)//&
	1094	" emis_tho_eq = ",emis_tho_eq
	1095
	1096	if (is_master)write(,)trim(rname)//&
	1097	"Latitude for diff albedo tholins"
	1098	tholateq=25. ! default value
	1099	call getin_p("tholateq",tholateq)
	1100	if (is_master)write(,)trim(rname)//&
	1101	" tholateq = ",tholateq
[3197]	1102	tholatn=-1.
	1103	tholats=0.
	1104	tholone=0.
	1105	tholonw=-1.
[3195]	1106	alb_tho_spe=0.1 ! default value
	1107	emis_tho_spe=1. ! default value
	1108	call getin_p(" tholatn",tholatn)
	1109	call getin_p(" tholats",tholats)
	1110	call getin_p(" tholone",tholone)
	1111	call getin_p(" tholonw",tholonw)
	1112	if (is_master)write(,)trim(rname)//&
	1113	" Values for special tholins albedo = ",tholatn,&
	1114	tholats,tholone,tholonw,alb_tho_spe,emis_tho_spe
	1115
	1116	if (is_master)write(,)trim(rname)//&
	1117	"Specific albedo"
	1118	spelon1=-180. ! default value
	1119	spelon2=180. ! default value
	1120	spelat1=-90. ! default value
	1121	spelat2=90. ! default value
	1122	specalb=.false. ! default value
	1123	if (is_master)write(,)trim(rname)//&
	1124	"albedo/emis spe "
	1125	albspe=0.1 ! default value
	1126	emispe=1. ! default value
	1127	call getin_p("spelon1",spelon1)
	1128	call getin_p("spelon2",spelon2)
	1129	call getin_p("spelat1",spelat1)
	1130	call getin_p("spelat2",spelat2)
	1131	call getin_p("specalb",specalb)
	1132	call getin_p("albspe",albspe)
	1133	call getin_p("emispe",emispe)
	1134
	1135	if (is_master)write(,)trim(rname)//&
	1136	" specific = ",specalb
	1137
	1138	!******************** TI ***************************
	1139
	1140	if (is_master)write(,)trim(rname)//&
	1141	"Change TI with time"
	1142	changeti=.false. ! default value
	1143	call getin_p("changeti",changeti)
	1144	if (is_master)write(,)trim(rname)//&
	1145	" changeti = ",changeti
	1146	changetid=.false. ! default value for diurnal TI
	1147	call getin_p("changetid",changetid)
	1148	if (is_master)write(,)trim(rname)//&
	1149	" changetid = ",changetid
	1150
	1151	if (is_master)write(,)trim(rname)//&
	1152	"IT N2 "
	1153	ITN2=800. ! default value
	1154	call getin_p("ITN2",ITN2)
	1155	if (is_master)write(,)trim(rname)//&
	1156	" ITN2 = ",ITN2
	1157	ITN2d=20. ! default value
	1158	call getin_p("ITN2d",ITN2d)
	1159	if (is_master)write(,)trim(rname)//&
	1160	" ITN2d = ",ITN2d
	1161
	1162	if (is_master)write(,)trim(rname)//&
	1163	"IT CH4"
	1164	ITCH4=800. ! default value
	1165	call getin_p("ITCH4",ITCH4)
	1166	if (is_master)write(,)trim(rname)//&
	1167	" ITCH4 = ",ITCH4
	1168	ITCH4d=20. ! default value
	1169	call getin_p("ITCH4d",ITCH4d)
	1170	if (is_master)write(,)trim(rname)//&
	1171	" ITCH4d = ",ITCH4d
	1172
	1173	if (is_master)write(,)trim(rname)//&
	1174	"IT H2O"
	1175	ITH2O=800. ! default value
	1176	call getin_p("ITH2O",ITH2O)
	1177	if (is_master)write(,)trim(rname)//&
	1178	" ITH2O = ",ITH2O
	1179	ITH2Od=20. ! default value
	1180	call getin_p("ITH2Od",ITH2Od)
	1181	if (is_master)write(,)trim(rname)//&
	1182	" ITH2Od = ",ITH2Od
	1183
[3372]	1184	!************ COOLING *************
	1185
	1186	alpha_top=5e-11 ! default value
	1187	call getin_p("alpha_top",alpha_top)
	1188	if (is_master)write(,)trim(rname)//&
	1189	" alpha_top = ",alpha_top
	1190	pref=0.02 ! default value
	1191	call getin_p("pref",pref)
	1192	if (is_master)write(,)trim(rname)//&
	1193	" pref = ",pref
	1194	deltap=0.1 ! default value
	1195	call getin_p("deltap",deltap)
	1196	if (is_master)write(,)trim(rname)//&
	1197	" deltap = ",deltap
	1198
[3184]	1199	!=================================
	1200	! TWOLAY scheme and NH3 cloudare left for retrocompatibility only,
	1201	! You should now use N-LAYER scheme (see below).
	1202
	1203	if (is_master) write(,)trim(rname)//&
	1204	": Radiatively active two-layer aerosols?"
	1205	aeroback2lay=.false. ! default value
	1206	call getin_p("aeroback2lay",aeroback2lay)
	1207	if (is_master) write(,)trim(rname)//": aeroback2lay = ",aeroback2lay
	1208
	1209	if (aeroback2lay.and.is_master) then
	1210	print*,'Warning : The TWOLAY AEROSOL scheme is deprecated and buggy...'
	1211	print*,'You should use the generic n-layer scheme (see aeronlay).'
	1212	endif
	1213
	1214	if (is_master) write(,)trim(rname)//": Radiatively active ammonia cloud?"
	1215	aeronh3=.false. ! default value
	1216	call getin_p("aeronh3",aeronh3)
	1217	if (is_master) write(,)trim(rname)//": aeronh3 = ",aeronh3
	1218
	1219	if (aeronh3.and.is_master) then
	1220	print*,'Warning : You are using specific NH3 cloud scheme ...'
	1221	print*,'You should use the generic n-layer scheme (see aeronlay).'
	1222	endif
	1223
	1224	if (is_master) write(,)trim(rname)//&
	1225	": TWOLAY AEROSOL: total optical depth "//&
	1226	"in the tropospheric layer (visible)"
	1227	obs_tau_col_tropo=8.D0
	1228	call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
	1229	if (is_master) write(,)trim(rname)//&
	1230	": obs_tau_col_tropo = ",obs_tau_col_tropo
	1231
	1232	if (is_master) write(,)trim(rname)//&
	1233	": TWOLAY AEROSOL: total optical depth "//&
	1234	"in the stratospheric layer (visible)"
	1235	obs_tau_col_strato=0.08D0
	1236	call getin_p("obs_tau_col_strato",obs_tau_col_strato)
	1237	if (is_master) write(,)trim(rname)//&
	1238	": obs_tau_col_strato = ",obs_tau_col_strato
	1239
	1240	if (is_master) write(,)trim(rname)//&
	1241	": TWOLAY AEROSOL: optprop_back2lay_vis?"
	1242	optprop_back2lay_vis = 'optprop_saturn_vis_n20.dat'
	1243	call getin_p("optprop_back2lay_vis",optprop_back2lay_vis)
	1244	if (is_master) write(,)trim(rname)//&
	1245	": optprop_back2lay_vis = ",trim(optprop_back2lay_vis)
	1246
	1247	if (is_master) write(,)trim(rname)//&
	1248	": TWOLAY AEROSOL: optprop_back2lay_ir?"
	1249	optprop_back2lay_ir = 'optprop_saturn_ir_n20.dat'
	1250	call getin_p("optprop_back2lay_ir",optprop_back2lay_ir)
	1251	if (is_master) write(,)trim(rname)//&
	1252	": optprop_back2lay_ir = ",trim(optprop_back2lay_ir)
[3197]	1253
[3184]	1254	if (is_master) write(,)trim(rname)//&
	1255	": TWOLAY AEROSOL: pres_bottom_tropo? in pa"
	1256	pres_bottom_tropo=66000.0
	1257	call getin_p("pres_bottom_tropo",pres_bottom_tropo)
	1258	if (is_master) write(,)trim(rname)//&
	1259	": pres_bottom_tropo = ",pres_bottom_tropo
	1260
	1261	if (is_master) write(,)trim(rname)//&
	1262	": TWOLAY AEROSOL: pres_top_tropo? in pa"
	1263	pres_top_tropo=18000.0
	1264	call getin_p("pres_top_tropo",pres_top_tropo)
	1265	if (is_master) write(,)trim(rname)//&
	1266	": pres_top_tropo = ",pres_top_tropo
	1267
	1268	if (is_master) write(,)trim(rname)//&
	1269	": TWOLAY AEROSOL: pres_bottom_strato? in pa"
	1270	pres_bottom_strato=2000.0
	1271	call getin_p("pres_bottom_strato",pres_bottom_strato)
	1272	if (is_master) write(,)trim(rname)//&
	1273	": pres_bottom_strato = ",pres_bottom_strato
	1274
	1275	! Sanity check
	1276	if (pres_bottom_strato .gt. pres_top_tropo) then
	1277	if(is_master) then
	1278	print*,'Error : TWOLAY AEROSOL, Please ensure that in callphys.def'
	1279	print*,'you have pres_top_tropo > pres_bottom_strato !'
	1280	endif
	1281	call abort_physic(rname," pres_top_tropo > pres_bottom_strato!",1)
	1282	endif
	1283
	1284	if (is_master) write(,)trim(rname)//&
	1285	": TWOLAY AEROSOL: pres_top_strato? in pa"
	1286	pres_top_strato=100.0
	1287	call getin_p("pres_top_strato",pres_top_strato)
	1288	if (is_master) write(,)trim(rname)//&
	1289	": pres_top_strato = ",pres_top_strato
	1290
	1291	if (is_master) write(,)trim(rname)//&
	1292	": TWOLAY AEROSOL: particle size in the ", &
	1293	"tropospheric layer, in meters"
	1294	size_tropo=2.e-6
	1295	call getin_p("size_tropo",size_tropo)
	1296	if (is_master) write(,)trim(rname)//": size_tropo = ",size_tropo
	1297
	1298	if (is_master) write(,)trim(rname)//&
	1299	": TWOLAY AEROSOL: particle size in the ", &
	1300	"stratospheric layer, in meters"
	1301	size_strato=1.e-7
	1302	call getin_p("size_strato",size_strato)
	1303	if (is_master) write(,)trim(rname)//": size_strato = ",size_strato
	1304
	1305	if (is_master) write(,)trim(rname)//&
	1306	": NH3 (thin) cloud: total optical depth"
	1307	tau_nh3_cloud=7.D0
	1308	call getin_p("tau_nh3_cloud",tau_nh3_cloud)
	1309	if (is_master) write(,)trim(rname)//": tau_nh3_cloud = ",tau_nh3_cloud
	1310
	1311	if (is_master) write(,)trim(rname)//": NH3 (thin) cloud pressure level"
	1312	pres_nh3_cloud=7.D0
	1313	call getin_p("pres_nh3_cloud",pres_nh3_cloud)
	1314	if (is_master) write(,)trim(rname)//": pres_nh3_cloud = ",pres_nh3_cloud
	1315
	1316	if (is_master) write(,)trim(rname)//": NH3 (thin) cloud: particle sizes"
	1317	size_nh3_cloud=3.e-6
	1318	call getin_p("size_nh3_cloud",size_nh3_cloud)
	1319	if (is_master) write(,)trim(rname)//": size_nh3_cloud = ",size_nh3_cloud
	1320
	1321	!=================================
	1322	! Generic N-LAYER aerosol scheme
	1323
	1324	if (is_master) write(,)trim(rname)//&
	1325	": Radiatively active generic n-layer aerosols?"
	1326	aeronlay=.false. ! default value
	1327	call getin_p("aeronlay",aeronlay)
	1328	if (is_master) write(,)trim(rname)//": aeronlay = ",aeronlay
	1329
	1330	if (is_master) write(,)trim(rname)//&
	1331	": Number of generic aerosols layers?"
	1332	nlayaero=1 ! default value
	1333	call getin_p("nlayaero",nlayaero)
	1334	! Avoid to allocate arrays of size 0
	1335	if (aeronlay .and. nlayaero.lt.1) then
	1336	if (is_master) then
	1337	print*, " You are trying to set no generic aerosols..."
	1338	print*, " Set aeronlay=.false. instead ! I abort."
	1339	endif
	1340	call abort_physic(rname,"no generic aerosols...",1)
	1341	endif
	1342	if (.not. aeronlay) nlayaero=1
	1343	if (is_master) write(,)trim(rname)//": nlayaero = ",nlayaero
	1344
	1345	! This is necessary, we just set the number of aerosol layers
[3197]	1346	IF(.NOT.ALLOCATED(aeronlay_tauref)) ALLOCATE(aeronlay_tauref(nlayaero))
	1347	IF(.NOT.ALLOCATED(aeronlay_lamref)) ALLOCATE(aeronlay_lamref(nlayaero))
	1348	IF(.NOT.ALLOCATED(aeronlay_choice)) ALLOCATE(aeronlay_choice(nlayaero))
	1349	IF(.NOT.ALLOCATED(aeronlay_pbot)) ALLOCATE(aeronlay_pbot(nlayaero))
	1350	IF(.NOT.ALLOCATED(aeronlay_ptop)) ALLOCATE(aeronlay_ptop(nlayaero))
	1351	IF(.NOT.ALLOCATED(aeronlay_sclhght)) ALLOCATE(aeronlay_sclhght(nlayaero))
	1352	IF(.NOT.ALLOCATED(aeronlay_size)) ALLOCATE(aeronlay_size(nlayaero))
	1353	IF(.NOT.ALLOCATED(aeronlay_nueff)) ALLOCATE(aeronlay_nueff(nlayaero))
	1354	IF(.NOT.ALLOCATED(optprop_aeronlay_ir)) ALLOCATE(optprop_aeronlay_ir(nlayaero))
	1355	IF(.NOT.ALLOCATED(optprop_aeronlay_vis)) ALLOCATE(optprop_aeronlay_vis(nlayaero))
[3184]	1356
	1357	if (is_master) write(,)trim(rname)//&
	1358	": Generic n-layer aerosols: Optical depth at reference wavelenght"
	1359	aeronlay_tauref=1.0E-1
	1360	call getin_p("aeronlay_tauref",aeronlay_tauref)
	1361	if (is_master) write(,)trim(rname)//&
	1362	": aeronlay_tauref = ",aeronlay_tauref
	1363
	1364	if (is_master) write(,)trim(rname)//&
	1365	": Generic n-layer aerosols: Reference wavelenght for optical depths (m)"
	1366	aeronlay_lamref=0.6E-6
	1367	call getin_p("aeronlay_lamref",aeronlay_lamref)
	1368	if (is_master) write(,)trim(rname)//&
	1369	": aeronlay_lamref = ",aeronlay_lamref
	1370
	1371	if (is_master) then
	1372	write(,)trim(rname)//&
	1373	": Generic n-layer aerosols: Vertical profile choice : "
	1374	write(,)trim(rname)//&
[3197]	1375	" (1) Tau btwn ptop and pbot follows atm. scale height"
[3184]	1376	write(,)trim(rname)//&
	1377	" or (2) Tau above pbot follows its own scale height"
	1378	endif
	1379	aeronlay_choice=1
	1380	call getin_p("aeronlay_choice",aeronlay_choice)
	1381	if (is_master) write(,)trim(rname)//&
	1382	": aeronlay_choice = ",aeronlay_choice
	1383
	1384	if (is_master) write(,)trim(rname)//&
	1385	": Generic n-layer aerosols: bottom pressures (Pa)"
	1386	aeronlay_pbot=2000.0
	1387	call getin_p("aeronlay_pbot",aeronlay_pbot)
	1388	if (is_master) write(,)trim(rname)//": aeronlay_pbot = ",aeronlay_pbot
	1389
	1390	if (is_master) write(,)trim(rname)//&
	1391	": Generic n-layer aerosols: (if choice=1) Top pressures (Pa) "
	1392	aeronlay_ptop=300000.0
	1393	call getin_p("aeronlay_ptop",aeronlay_ptop)
	1394	if (is_master) write(,)trim(rname)//": aeronlay_ptop = ",aeronlay_ptop
	1395
	1396	if (is_master) write(,)trim(rname)//&
	1397	": Generic n-layer aerosols: (if choice=2) Scale height / atm. scale height"
	1398	aeronlay_sclhght=0.2
	1399	call getin_p("aeronlay_sclhght",aeronlay_sclhght)
	1400	if (is_master) write(,)trim(rname)//&
	1401	": aeronlay_sclhght = ",aeronlay_sclhght
	1402
	1403	if (is_master) write(,)trim(rname)//&
	1404	": Generic n-layer aerosols: particles effective radii (m)"
	1405	aeronlay_size=1.e-6
	1406	call getin_p("aeronlay_size",aeronlay_size)
	1407	if (is_master) write(,)trim(rname)//": aeronlay_size = ",aeronlay_size
	1408
	1409	if (is_master) write(,)trim(rname)//&
	1410	": Generic n-layer aerosols: particles radii effective variance"
	1411	aeronlay_nueff=0.1
	1412	call getin_p("aeronlay_nueff",aeronlay_nueff)
	1413	if (is_master) write(,)trim(rname)//": aeronlay_nueff = ",aeronlay_nueff
	1414
	1415	if (is_master) write(,)trim(rname)//&
	1416	": Generic n-layer aerosols: VIS optical properties file"
	1417	optprop_aeronlay_vis = 'optprop_saturn_vis_n20.dat'
	1418	call getin_p("optprop_aeronlay_vis",optprop_aeronlay_vis)
	1419	if (is_master) write(,)trim(rname)//&
	1420	": optprop_aeronlay_vis = ",optprop_aeronlay_vis
	1421
	1422	if (is_master) write(,)trim(rname)//&
	1423	": Generic n-layer aerosols: IR optical properties file"
	1424	optprop_aeronlay_ir = 'optprop_saturn_ir_n20.dat'
	1425	call getin_p("optprop_aeronlay_ir",optprop_aeronlay_ir)
	1426	if (is_master) write(,)trim(rname)//&
	1427	": optprop_aeronlay_ir = ",optprop_aeronlay_ir
	1428
[3197]	1429
[3184]	1430	!=================================
	1431
	1432	if (is_master) write(,)trim(rname)//&
	1433	": Is the variable gas species radiatively active?"
	1434	Tstrat=167.0
	1435	varactive=.false.
	1436	call getin_p("varactive",varactive)
	1437	if (is_master) write(,)trim(rname)//": varactive = ",varactive
	1438
	1439	if (is_master) write(,)trim(rname)//&
	1440	": Is the variable gas species distribution set?"
	1441	varfixed=.false.
	1442	call getin_p("varfixed",varfixed)
	1443	if (is_master) write(,)trim(rname)//": varfixed = ",varfixed
	1444
	1445	if (is_master) write(,)trim(rname)//&
	1446	": What is the saturation % of the variable species?"
	1447	satval=0.8
	1448	call getin_p("satval",satval)
	1449	if (is_master) write(,)trim(rname)//": satval = ",satval
	1450
	1451
	1452	! Test of incompatibility:
	1453	! if varactive, then varfixed should be false
	1454	if (varactive.and.varfixed) then
	1455	call abort_physic(rname,'if varactive, varfixed must be OFF!',1)
	1456	endif
	1457
	1458	if (is_master) write(,)trim(rname)//": Gravitationnal sedimentation ?"
	1459	sedimentation=.false. ! default value
	1460	call getin_p("sedimentation",sedimentation)
	1461	if (is_master) write(,)trim(rname)//": sedimentation = ",sedimentation
	1462
	1463	if (is_master) write(,)trim(rname)//": Generic Condensation of tracers ?"
[3197]	1464	generic_condensation=.false. !default value
[3184]	1465	call getin_p("generic_condensation",generic_condensation)
	1466	if (is_master) write(,)trim(rname)//": generic_condensation = ",generic_condensation
[3197]	1467
[3184]	1468	if (is_master) write(,)trim(rname)//": Spectral Dependant albedo ?"
	1469	albedo_spectral_mode=.false. ! default value
	1470	call getin_p("albedo_spectral_mode",albedo_spectral_mode)
	1471	if (is_master) write(,)trim(rname)//&
	1472	": albedo_spectral_mode = ",albedo_spectral_mode
	1473
	1474	if (is_master) then
	1475	write(,)trim(rname)//": Snow albedo ?"
	1476	write(,)trim(rname)//": If albedo_spectral_mode=.true., then this "
	1477	write(,)trim(rname)//": corresponds to the 0.5 microns snow albedo."
	1478	endif
	1479	albedosnow=0.5 ! default value
	1480	call getin_p("albedosnow",albedosnow)
	1481	if (is_master) write(,)trim(rname)//": albedosnow = ",albedosnow
[3197]	1482
[3184]	1483	if (is_master) write(,)trim(rname)//": N2 ice albedo ?"
	1484	albedon2ice=0.5 ! default value
	1485	call getin_p("albedon2ice",albedon2ice)
	1486	if (is_master) write(,)trim(rname)//": albedon2ice = ",albedon2ice
	1487
	1488	if (is_master) write(,)trim(rname)//": Maximum ice thickness ?"
	1489	maxicethick=2.0 ! default value
	1490	call getin_p("maxicethick",maxicethick)
	1491	if (is_master) write(,)trim(rname)//": maxicethick = ",maxicethick
	1492
	1493	if (is_master) write(,)trim(rname)//": Minimum eddy mix coeff in 1D ?"
	1494	kmixmin=1.0e-2 ! default value
	1495	call getin_p("kmixmin",kmixmin)
	1496	if (is_master) write(,)trim(rname)//": kmixmin = ",kmixmin
	1497
	1498	if (is_master) write(,)'Predefined Cp from dynamics is ',cpp,'J kg^-1 K^-1'
	1499	if (is_master) write(,)'Predefined Mg from dynamics is ',mugaz,'amu'
	1500
	1501	force_cpp=.false. ! default value
	1502	call getin_p("force_cpp",force_cpp)
	1503	if (force_cpp) then
	1504	if (is_master) write(,)trim(rname)//": force_cpp = ",force_cpp
	1505	if (is_master) write(,)trim(rname)//": force_cpp is deprecated.",&
	1506	"Set cpp_mugaz_mode=1 in callfis to emulate force_cpp=.true."
	1507	call abort_physic(rname,"Anyway, you need to set force_cpp=.false. to continue.",1)
	1508	endif
	1509
	1510	if (is_master) write(,)trim(rname)//&
	1511	": where do you want your cpp/mugaz value to come from?",&
	1512	"=> 0: dynamics (3d), 1: forced in callfis (1d), 2: computed from gases.def (1d)?"
	1513	cpp_mugaz_mode = 0 ! default value
	1514	call getin_p("cpp_mugaz_mode",cpp_mugaz_mode)
	1515	if (is_master) write(,)trim(rname)//": cpp_mugaz_mode = ",cpp_mugaz_mode
	1516
[3422]	1517
	1518	! Test of incompatibility:
	1519
[3439]	1520	if ((.not.tracer).and.(haze)) then
	1521	call abort_physic(rname, 'if haze are on, tracers must be on!', 1)
	1522	endif
[3477]	1523	if (haze_proffix.and.sedimentation) then
	1524	call abort_physic(rname, 'if haze profile is set, sedimentation must be deactivated', 1)
	1525	endif
[3455]	1526	if (callmolvis.and..not.callconduct) then
	1527	call abort_physic(rname, 'if callmolvis is set, callconduct must be true', 1)
	1528	endif
[3474]	1529	if (glaflow.and..not.fast) then
	1530	call abort_physic(rname, 'if glaflow is set, fast must be true', 1)
	1531	endif
[3439]	1532	if (paleo.and..not.fast) then
	1533	call abort_physic(rname, 'if paleo is set, fast must be true', 1)
[3455]	1534	endif
[3439]	1535	if ((haze_proffix.or.haze_radproffix).and..not.aerohaze) then
	1536	call abort_physic(rname, 'for now, haze/rad proffix only works w aerohaze=T', 1)
[3455]	1537	endif
[3539]	1538	if (carbox.and.condcosurf.and.no_n2frost) then
[3439]	1539	call abort_physic(rname, "CO surface condensation and no_n2frost are both active which may not be relevant", 1)
	1540	end if
	1541
[3184]	1542	if ((cpp_mugaz_mode >= 1).and.(is_master).and.(ngrid>1)) then
	1543	write(,)' !!!!! Be aware that having different values for cpp and mugaz in the dynamics and physics'
	1544	write(,)' in 3D can result in very pathological behavior. You have been warned !!!!!'
	1545	else if ((cpp_mugaz_mode == 0).and.(is_master).and.(ngrid==1)) then
	1546	! for this specific 1D error we will remove run.def before aborting JL22
	1547	call system("rm -rf run.def")
	1548	call abort_physic(rname,"cpp_mugaz_mode must be >= 1 in 1d",1)
	1549	endif
	1550
	1551	if (cpp_mugaz_mode == 1) then
	1552	mugaz = -99999.
	1553	if (is_master) write(,)trim(rname)//&
	1554	": MEAN MOLECULAR MASS in g mol-1 ?"
	1555	call getin_p("mugaz",mugaz)
	1556	IF (mugaz.eq.-99999.) THEN
	1557	call abort_physic(rname,"mugaz must be set if cpp_mugaz_mode = 1",1)
	1558	ENDIF
	1559	cpp = -99999.
	1560	if (is_master) write(,)trim(rname)//&
	1561	": SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
	1562	call getin_p("cpp",cpp)
	1563	IF (cpp.eq.-99999.) THEN
	1564	PRINT *, "cpp must be set if cpp_mugaz_mode = 1"
	1565	STOP
	1566	ENDIF
	1567	if (is_master) write(,)'New Cp from callfis is ',cpp,'J kg^-1 K^-1'
	1568	if (is_master) write(,)'New Mg from callfis is ',mugaz,'amu'
[3197]	1569
[3184]	1570	endif ! of if (cpp_mugaz_mode == 1)
	1571	call su_gases
	1572	call calc_cpp_mugaz
	1573
	1574	if (is_master) then
	1575	PRINT*,'--------------------------------------------'
	1576	PRINT*
	1577	PRINT*
	1578	endif
	1579	ELSE
	1580	call abort_physic(rname,'Cannot read file callphys.def. Is it here ?',1)
	1581	ENDIF ! of IF(iscallphys)
	1582
	1583	if (is_master) then
	1584	PRINT*
	1585	PRINT*,'inifis: daysec',daysec
	1586	PRINT*
	1587	PRINT*,'inifis: The radiative transfer is computed:'
	1588	PRINT*,' each ',iradia,' physical time-step'
	1589	PRINT,' or each ',iradiadtphys,' seconds'
	1590	PRINT*
	1591	endif
	1592
	1593	!-----------------------------------------------------------------------
	1594	! Some more initialization:
	1595	! ------------------------
	1596
	1597	! Initializations for comgeomfi_h
	1598	#ifndef MESOSCALE
	1599	totarea=SSUM(ngrid,parea,1)
	1600	call planetwide_sumval(parea,totarea_planet)
	1601
	1602	!! those are defined in comdiurn_h.F90
	1603	IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
	1604	IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
	1605	IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
	1606	IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
	1607
	1608	DO ig=1,ngrid
	1609	sinlat(ig)=sin(plat(ig))
	1610	coslat(ig)=cos(plat(ig))
	1611	sinlon(ig)=sin(plon(ig))
	1612	coslon(ig)=cos(plon(ig))
	1613	ENDDO
	1614	#endif
	1615	! initialize variables in radinc_h
	1616	call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
	1617
	1618	! initialize variables and allocate arrays in radcommon_h
	1619	call ini_radcommon_h(naerkind)
[3197]	1620
[3184]	1621	! allocate "comsoil_h" arrays
	1622	call ini_comsoil_h(ngrid)
[3197]	1623
[3184]	1624	END SUBROUTINE inifis
	1625
	1626	END MODULE inifis_mod

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

Download in other formats:

Original Format