Context Navigation

inifis_mod.F90 @ 3652

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

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

Download in other formats:

Original Format