Context Navigation

inifis_mod.F90 @ 3607

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

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

Download in other formats:

Original Format