Context Navigation

inifis_mod.F90 @ 3579

Last change on this file since 3579 was 3437, checked in by emillour, 9 months ago

Generic PCM Integrate update from newton branch by LT:

remove call to old newtrelax routine and replace it with a new one, newtcool_MOCHA, tuned to the MOCHA intercomparison project

File size: 49.7 KB

Rev	Line
[1524]	1	MODULE inifis_mod
	2	IMPLICIT NONE
[253]	3
[1524]	4	CONTAINS
[253]	5
[1524]	6	SUBROUTINE inifis(ngrid,nlayer,nq, &
[1525]	7	day_ini,pdaysec,nday,ptimestep, &
[1524]	8	plat,plon,parea, &
	9	prad,pg,pr,pcpp)
	10
[1682]	11	use init_print_control_mod, only: init_print_control
[1529]	12	use radinc_h, only: ini_radinc_h, naerkind
[2972]	13	use radcommon_h, only: ini_radcommon_h
[2954]	14	use radii_mod, only: radfixed, Nmix_co2
[1524]	15	use datafile_mod, only: datadir
	16	use comdiurn_h, only: sinlat, coslat, sinlon, coslon
[1542]	17	use comgeomfi_h, only: totarea, totarea_planet
[1538]	18	use comsoil_h, only: ini_comsoil_h, nsoilmx, lay1_soil, alpha_soil
[1525]	19	use time_phylmdz_mod, only: ecritphy,day_step,iphysiq, &
	20	init_time, daysec, dtphys
[1524]	21	use comcstfi_mod, only: rad, cpp, g, r, rcp, &
	22	mugaz, pi, avocado
	23	use planete_mod, only: nres
	24	use planetwide_mod, only: planetwide_sumval
	25	use callkeys_mod
[2958]	26	use wstats_mod, only: callstats
[2299]	27	use ioipsl_getin_p_mod, only : getin_p
[2583]	28	use mod_phys_lmdz_para, only : is_parallel, is_master, bcast
[3437]	29	use newton_cooling_hotJ, only: planetary_suffix
[1524]	30
[135]	31	!=======================================================================
	32	!
	33	! purpose:
	34	! -------
	35	!
	36	! Initialisation for the physical parametrisations of the LMD
[305]	37	! Generic Model.
[135]	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	! -------------
[1524]	66	use datafile_mod, only: datadir
	67	use ioipsl_getin_p_mod, only: getin_p
	68	IMPLICIT NONE
[1384]	69
[135]	70
	71
[1524]	72	REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep
[1525]	73	INTEGER,INTENT(IN) :: nday
[1524]	74	INTEGER,INTENT(IN) :: ngrid,nlayer,nq
	75	REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid)
	76	integer,intent(in) :: day_ini
[2583]	77
	78	INTEGER ig
	79	CHARACTER(len=20) :: rname="inifis" ! routine name, for messages
[135]	80
[1524]	81	EXTERNAL iniorbit,orbite
	82	EXTERNAL SSUM
	83	REAL SSUM
[135]	84
[1682]	85	! Initialize flags lunout, prt_level, debug (in print_control_mod)
	86	CALL init_print_control
	87
[1524]	88	! initialize constants in comcstfi_mod
	89	rad=prad
	90	cpp=pcpp
	91	g=pg
	92	r=pr
	93	rcp=r/cpp
[2078]	94	mugaz=8.314*1000./pr ! dummy init
[1524]	95	pi=2.*asin(1.)
	96	avocado = 6.02214179e23 ! added by RW
[135]	97
[1524]	98	! Initialize some "temporal and calendar" related variables
[1832]	99	#ifndef MESOSCALE
[1525]	100	CALL init_time(day_ini,pdaysec,nday,ptimestep)
[1832]	101	#endif
[135]	102
[1525]	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
[135]	109
[1669]	110	! do we read a startphy.nc file? (default: .true.)
	111	call getin_p("startphy_file",startphy_file)
	112
[135]	113	! --------------------------------------------------------------
	114	! Reading the "callphys.def" file controlling some key options
	115	! --------------------------------------------------------------
	116
[2583]	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
[1524]	124	IF(iscallphys) THEN
[135]	125
[2583]	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:"
[1524]	135	! default 'datadir' is set in "datadir_mod"
	136	call getin_p("datadir",datadir) ! default path
[2583]	137	if (is_master) write(,) trim(rname)//": datadir = ",trim(datadir)
[374]	138
[2583]	139	if (is_master) write(,) trim(rname)//&
	140	": Run with or without tracer transport ?"
[1524]	141	tracer=.false. ! default value
	142	call getin_p("tracer",tracer)
[2583]	143	if (is_master) write(,) trim(rname)//": tracer = ",tracer
[135]	144
[2583]	145	if (is_master) write(,) trim(rname)//&
	146	": Run with or without atm mass update "//&
	147	" due to tracer evaporation/condensation?"
[1524]	148	mass_redistrib=.false. ! default value
	149	call getin_p("mass_redistrib",mass_redistrib)
[2583]	150	if (is_master) write(,) trim(rname)//&
	151	": mass_redistrib = ",mass_redistrib
[728]	152
[2583]	153	if (is_master) then
	154	write(,) trim(rname)//": Diurnal cycle ?"
	155	write(,) trim(rname)//&
	156	": (if diurnal=false, diurnal averaged solar heating)"
	157	endif
[1524]	158	diurnal=.true. ! default value
	159	call getin_p("diurnal",diurnal)
[2583]	160	if (is_master) write(,) trim(rname)//": diurnal = ",diurnal
[135]	161
[2583]	162	if (is_master) then
	163	write(,) trim(rname)//": Seasonal cycle ?"
	164	write(,) trim(rname)//&
	165	": (if season=false, Ls stays constant, to value ", &
	166	"set in 'start'"
	167	endif
[1524]	168	season=.true. ! default value
	169	call getin_p("season",season)
[2583]	170	if (is_master) write(,) trim(rname)//": season = ",season
[2300]	171
[2583]	172	if (is_master) write(,) trim(rname)//&
	173	": No seasonal cycle: initial day to lock the run during restart"
[2300]	174	noseason_day=0.0 ! default value
	175	call getin_p("noseason_day",noseason_day)
[2583]	176	if (is_master) write(,) trim(rname)//": noseason_day=", noseason_day
[135]	177
[2583]	178	if (is_master) write(,) trim(rname)//": Tidally resonant rotation ?"
[1524]	179	tlocked=.false. ! default value
	180	call getin_p("tlocked",tlocked)
[2583]	181	if (is_master) write(,) trim(rname)//": tlocked = ",tlocked
[135]	182
[2583]	183	if (is_master) write(,) trim(rname)//": Saturn ring shadowing ?"
[1524]	184	rings_shadow = .false.
	185	call getin_p("rings_shadow", rings_shadow)
[2583]	186	if (is_master) write(,) trim(rname)//": rings_shadow = ", rings_shadow
[1133]	187
[2583]	188	if (is_master) write(,) trim(rname)//&
	189	": Compute latitude-dependent gravity field?"
[1524]	190	oblate = .false.
	191	call getin_p("oblate", oblate)
[2583]	192	if (is_master) write(,) trim(rname)//": oblate = ", oblate
[1194]	193
[2583]	194	if (is_master) write(,) trim(rname)//": Flattening of the planet (a-b)/a "
[1524]	195	flatten = 0.0
	196	call getin_p("flatten", flatten)
[2583]	197	if (is_master) write(,) trim(rname)//": flatten = ", flatten
[1174]	198
[1194]	199
[2583]	200	if (is_master) write(,) trim(rname)//": Needed if oblate=.true.: J2"
[1524]	201	J2 = 0.0
	202	call getin_p("J2", J2)
[2583]	203	if (is_master) write(,) trim(rname)//": J2 = ", J2
[1194]	204
[2583]	205	if (is_master) write(,) trim(rname)//&
	206	": Needed if oblate=.true.: Planet mass (*1e24 kg)"
[1524]	207	MassPlanet = 0.0
	208	call getin_p("MassPlanet", MassPlanet)
[2583]	209	if (is_master) write(,) trim(rname)//": MassPlanet = ", MassPlanet
[1194]	210
[2583]	211	if (is_master) write(,) trim(rname)//&
	212	": Needed if oblate=.true.: Planet mean radius (m)"
[1524]	213	Rmean = 0.0
	214	call getin_p("Rmean", Rmean)
[2583]	215	if (is_master) write(,) trim(rname)//": Rmean = ", Rmean
[1194]	216
[135]	217	! Test of incompatibility:
	218	! if tlocked, then diurnal should be false
[1524]	219	if (tlocked.and.diurnal) then
[2583]	220	call abort_physic(rname,"If diurnal=true, we should turn off tlocked",1)
[1524]	221	endif
[135]	222
[2583]	223	if (is_master) write(,) trim(rname)//": Tidal resonance ratio ?"
[1524]	224	nres=0 ! default value
	225	call getin_p("nres",nres)
[2583]	226	if (is_master) write(,) trim(rname)//": nres = ",nres
[135]	227
[2583]	228	if (is_master) write(,) trim(rname)//&
	229	": Write some extra output to the screen ?"
[1524]	230	lwrite=.false. ! default value
	231	call getin_p("lwrite",lwrite)
[2583]	232	if (is_master) write(,) trim(rname)//": lwrite = ",lwrite
[135]	233
[2583]	234	if (is_master) write(,) trim(rname)//&
	235	": Save statistics in file stats.nc ?"
[1524]	236	callstats=.true. ! default value
	237	call getin_p("callstats",callstats)
[2583]	238	if (is_master) write(,) trim(rname)//": callstats = ",callstats
[135]	239
[2583]	240	if (is_master) write(,) trim(rname)//&
	241	": Test energy conservation of model physics ?"
[1524]	242	enertest=.false. ! default value
	243	call getin_p("enertest",enertest)
[2583]	244	if (is_master) write(,) trim(rname)//": enertest = ",enertest
[135]	245
[2583]	246	if (is_master) write(,) trim(rname)//": call radiative transfer ?"
[1524]	247	callrad=.true. ! default value
	248	call getin_p("callrad",callrad)
[2583]	249	if (is_master) write(,) trim(rname)//": callrad = ",callrad
[135]	250
[2583]	251	if (is_master) write(,) trim(rname)//&
	252	": call correlated-k radiative transfer ?"
[1524]	253	corrk=.true. ! default value
	254	call getin_p("corrk",corrk)
[2583]	255	if (is_master) write(,) trim(rname)//": corrk = ",corrk
[135]	256
[2583]	257	if (is_master) write(,) trim(rname)//&
	258	": prohibit calculations outside corrk T grid?"
[1524]	259	strictboundcorrk=.true. ! default value
	260	call getin_p("strictboundcorrk",strictboundcorrk)
[2583]	261	if (is_master) write(,) trim(rname)//&
	262	": strictboundcorrk = ",strictboundcorrk
[3235]	263
	264	if (is_master) write(,) trim(rname)//&
	265	": call thermosphere ?"
	266	callthermos=.false. ! default value
	267	call getin_p("callthermos",callthermos)
	268	if (is_master) write(,) trim(rname)//&
	269	": callthermos = ",callthermos
	270	if(callthermos) then
	271	if (is_master) write(,) trim(rname)//&
	272	": flux from thermosphere ? W/m2"
	273	phitop = 0.0 ! default value
	274	call getin_p("phitop",phitop)
	275	if (is_master) write(,) trim(rname)//&
	276	": phitop = ",phitop
	277	if (is_master) write(,) trim(rname)//&
	278	": Thickness of conduction ? m"
	279	zztop = 1000.
	280	call getin_p("zztop",zztop)
	281	if (is_master) write(,) trim(rname)//&
	282	": zztop = ",zztop
	283	if (is_master) write(,) trim(rname)//&
	284	": Force conduction ? "
	285	force_conduction=.false. ! default value
	286	call getin_p("force_conduction",force_conduction)
	287	if (is_master) write(,) trim(rname)//&
	288	": force_conduction = ",force_conduction
	289	if(force_conduction) then
	290	if (is_master) write(,) trim(rname)//&
	291	": a_coeff ?"
	292	a_coeff = 0.0 ! default value
	293	call getin_p("a_coeff",a_coeff)
	294	if (is_master) write(,) trim(rname)//&
	295	": a_coeff = ",a_coeff
	296	if (is_master) write(,) trim(rname)//&
	297	": s_coeff ?"
	298	s_coeff = 0.0 ! default value
	299	call getin_p("s_coeff",s_coeff)
	300	if (is_master) write(,) trim(rname)//&
	301	": s_coeff = ",s_coeff
	302	endif
	303	endif
[1145]	304
[2583]	305	if (is_master) write(,) trim(rname)//&
[2655]	306	": prohibit calculations outside CIA T grid?"
	307	strictboundcia=.true. ! default value
	308	call getin_p("strictboundcia",strictboundcia)
[2633]	309	if (is_master) write(,) trim(rname)//&
[2655]	310	": strictboundcia = ",strictboundcia
[2633]	311
	312	if (is_master) write(,) trim(rname)//&
[2583]	313	": Minimum atmospheric temperature for Planck function integration ?"
[2283]	314	tplanckmin=30.0 ! default value
	315	call getin_p("tplanckmin",tplanckmin)
[2583]	316	if (is_master) write(,) trim(rname)//": tplanckmin = ",tplanckmin
[2283]	317
[2583]	318	if (is_master) write(,) trim(rname)//&
	319	": Maximum atmospheric temperature for Planck function integration ?"
[2283]	320	tplanckmax=1500.0 ! default value
	321	call getin_p("tplanckmax",tplanckmax)
[2583]	322	if (is_master) write(,) trim(rname)//": tplanckmax = ",tplanckmax
[2283]	323
[2583]	324	if (is_master) write(,) trim(rname)//&
	325	": Temperature step for Planck function integration ?"
[2283]	326	dtplanck=0.1 ! default value
	327	call getin_p("dtplanck",dtplanck)
[2583]	328	if (is_master) write(,) trim(rname)//": dtplanck = ",dtplanck
[2283]	329
[2583]	330	if (is_master) write(,) trim(rname)//&
	331	": call gaseous absorption in the visible bands?"//&
	332	" (matters only if callrad=T)"
[1524]	333	callgasvis=.false. ! default value
	334	call getin_p("callgasvis",callgasvis)
[2583]	335	if (is_master) write(,) trim(rname)//": callgasvis = ",callgasvis
[538]	336
[2583]	337	if (is_master) write(,) trim(rname)//&
	338	": call continuum opacities in radiative transfer ?"//&
	339	" (matters only if callrad=T)"
[1524]	340	continuum=.true. ! default value
	341	call getin_p("continuum",continuum)
[2583]	342	if (is_master) write(,) trim(rname)//": continuum = ",continuum
[538]	343
[2583]	344	if (is_master) write(,) trim(rname)//": version for H2H2 CIA file ?"
[2245]	345	versH2H2cia=2011 ! default value (should be 2018 but retrocompatibility first)
	346	call getin_p("versH2H2cia",versH2H2cia)
[2583]	347	if (is_master) write(,) trim(rname)//": versH2H2cia = ",versH2H2cia
[2245]	348	! Sanity check
	349	if (versH2H2cia.ne.2011 .and. versH2H2cia.ne.2018) then
[2583]	350	call abort_physic(rname,"Error: Choose a correct value (2011 or 2018) for versH2H2cia !",1)
[2655]	351	endif
	352
	353	if (is_master) write(,) trim(rname)//&
	354	": CIA - normal or equilibrium ortho-para mixture for H2?"
	355	H2orthopara_mixture = 'normal'
	356	call getin_p("H2orthopara_mixture",H2orthopara_mixture)
	357	if (is_master) write(,)trim(rname)//&
	358	": H2orthopara_mixture = ",trim(H2orthopara_mixture)
[2245]	359
[2583]	360	if (is_master) write(,) trim(rname)//&
	361	": call turbulent vertical diffusion ?"
[1524]	362	calldifv=.true. ! default value
	363	call getin_p("calldifv",calldifv)
[2583]	364	if (is_master) write(,) trim(rname)//": calldifv = ",calldifv
[135]	365
[2583]	366	if (is_master) write(,) trim(rname)//": use turbdiff instead of vdifc ?"
[1524]	367	UseTurbDiff=.true. ! default value
	368	call getin_p("UseTurbDiff",UseTurbDiff)
[2583]	369	if (is_master) write(,) trim(rname)//": UseTurbDiff = ",UseTurbDiff
[596]	370
[2583]	371	if (is_master) write(,) trim(rname)//": call convective adjustment ?"
[1524]	372	calladj=.true. ! default value
	373	call getin_p("calladj",calladj)
[2583]	374	if (is_master) write(,) trim(rname)//": calladj = ",calladj
[135]	375
[2583]	376	if (is_master) write(,)trim(rname)//&
	377	": call Lott's non-oro GWs parameterisation scheme ?"
[2299]	378	calllott_nonoro=.false. ! default value
	379	call getin_p("calllott_nonoro",calllott_nonoro)
[2583]	380	if (is_master) write(,)trim(rname)//&
	381	": calllott_nonoro = ",calllott_nonoro
[2299]	382
[2583]	383	if (is_master) write(,) trim(rname)//": call thermal plume model ?"
[2060]	384	calltherm=.false. ! default value
	385	call getin_p("calltherm",calltherm)
[2583]	386	if (is_master) write(,) trim(rname)//": calltherm = ",calltherm
[2060]	387
[2583]	388	if (is_master) write(,) trim(rname)//": call CO2 condensation ?"
[1524]	389	co2cond=.false. ! default value
	390	call getin_p("co2cond",co2cond)
[2583]	391	if (is_master) write(,) trim(rname)//": co2cond = ",co2cond
[650]	392	! Test of incompatibility
[1524]	393	if (co2cond.and.(.not.tracer)) then
[2583]	394	call abort_physic(rname,"Error: We need a CO2 ice tracer to condense CO2!",1)
[1524]	395	endif
[716]	396
[2583]	397	if (is_master) write(,) trim(rname)//": CO2 supersaturation level ?"
[1524]	398	co2supsat=1.0 ! default value
	399	call getin_p("co2supsat",co2supsat)
[2583]	400	if (is_master) write(,) trim(rname)//": co2supsat = ",co2supsat
[253]	401
[2583]	402	if (is_master) write(,) trim(rname)//&
	403	": Radiative timescale for Newtonian cooling (in Earth days)?"
[1524]	404	tau_relax=30. ! default value
	405	call getin_p("tau_relax",tau_relax)
[2583]	406	if (is_master) write(,) trim(rname)//": tau_relax = ",tau_relax
[1524]	407	tau_relax=tau_relax243600 ! convert Earth days --> seconds
[253]	408
[3437]	409	if (is_master) write(,) trim(rname)//&
	410	": planetary suffix for loading file in Newtonian cooling HotJ ?"
	411	planetary_suffix='default_value' ! default value
	412	call getin_p("planetary_suffix",planetary_suffix)
	413	if (is_master) write(,) trim(rname)//": planetary_suffix",planetary_suffix
	414
[2583]	415	if (is_master) write(,)trim(rname)//&
	416	": call thermal conduction in the soil ?"
[1524]	417	callsoil=.true. ! default value
	418	call getin_p("callsoil",callsoil)
[2583]	419	if (is_master) write(,) trim(rname)//": callsoil = ",callsoil
[135]	420
[2583]	421	if (is_master) write(,)trim(rname)//&
	422	": Rad transfer is computed every iradia", &
	423	" physical timestep"
[1524]	424	iradia=1 ! default value
	425	call getin_p("iradia",iradia)
[2583]	426	if (is_master) write(,)trim(rname)//": iradia = ",iradia
[135]	427
[2583]	428	if (is_master) write(,)trim(rname)//": Rayleigh scattering ?"
[1524]	429	rayleigh=.false.
	430	call getin_p("rayleigh",rayleigh)
[2583]	431	if (is_master) write(,)trim(rname)//": rayleigh = ",rayleigh
[135]	432
[2583]	433	if (is_master) write(,) trim(rname)//&
	434	": Use blackbody for stellar spectrum ?"
[1524]	435	stelbbody=.false. ! default value
	436	call getin_p("stelbbody",stelbbody)
[2583]	437	if (is_master) write(,) trim(rname)//": stelbbody = ",stelbbody
[135]	438
[2583]	439	if (is_master) write(,) trim(rname)//": Stellar blackbody temperature ?"
[1524]	440	stelTbb=5800.0 ! default value
	441	call getin_p("stelTbb",stelTbb)
[2583]	442	if (is_master) write(,) trim(rname)//": stelTbb = ",stelTbb
[253]	443
[2583]	444	if (is_master) write(,) trim(rname)//": Output mean OLR in 1D?"
[1524]	445	meanOLR=.false.
	446	call getin_p("meanOLR",meanOLR)
[2583]	447	if (is_master) write(,) trim(rname)//": meanOLR = ",meanOLR
[135]	448
[2583]	449	if (is_master) write(,)trim(rname)//": Output spectral OLR in 3D?"
[1524]	450	specOLR=.false.
	451	call getin_p("specOLR",specOLR)
[2583]	452	if (is_master) write(,)trim(rname)//": specOLR = ",specOLR
[135]	453
[2583]	454	if (is_master) write(,)trim(rname)//&
	455	": Output diagnostic optical thickness attenuation exp(-tau)?"
[2131]	456	diagdtau=.false.
	457	call getin_p("diagdtau",diagdtau)
[2583]	458	if (is_master) write(,)trim(rname)//": diagdtau = ",diagdtau
[2131]	459
[2583]	460	if (is_master) write(,)trim(rname)//": Operate in kastprof mode?"
[1524]	461	kastprof=.false.
	462	call getin_p("kastprof",kastprof)
[2583]	463	if (is_master) write(,)trim(rname)//": kastprof = ",kastprof
[253]	464
[2583]	465	if (is_master) write(,)trim(rname)//&
	466	": Uniform absorption in radiative transfer?"
[1524]	467	graybody=.false.
	468	call getin_p("graybody",graybody)
[2583]	469	if (is_master) write(,)trim(rname)//": graybody = ",graybody
[253]	470
[1538]	471	! Soil model
[2583]	472	if (is_master) write(,)trim(rname)//&
	473	": Number of sub-surface layers for soil scheme?"
[1538]	474	! default value of nsoilmx set in comsoil_h
	475	call getin_p("nsoilmx",nsoilmx)
[2583]	476	if (is_master) write(,)trim(rname)//": nsoilmx=",nsoilmx
[1538]	477
[2583]	478	if (is_master) write(,)trim(rname)//&
	479	": Thickness of topmost soil layer (m)?"
[1538]	480	! default value of lay1_soil set in comsoil_h
	481	call getin_p("lay1_soil",lay1_soil)
[2583]	482	if (is_master) write(,)trim(rname)//": lay1_soil = ",lay1_soil
[1538]	483
[2583]	484	if (is_master) write(,)trim(rname)//&
	485	": Coefficient for soil layer thickness distribution?"
[1538]	486	! default value of alpha_soil set in comsoil_h
	487	call getin_p("alpha_soil",alpha_soil)
[2583]	488	if (is_master) write(,)trim(rname)//": alpha_soil = ",alpha_soil
[1538]	489
[1297]	490	! Slab Ocean
[2583]	491	if (is_master) write(,) trim(rname)//": Use slab-ocean ?"
[1524]	492	ok_slab_ocean=.false. ! default value
	493	call getin_p("ok_slab_ocean",ok_slab_ocean)
[2583]	494	if (is_master) write(,) trim(rname)//": ok_slab_ocean = ",ok_slab_ocean
[1529]	495	! Sanity check: for now slab oncean only works in serial mode
[3100]	496	! if (ok_slab_ocean.and.is_parallel) then
	497	! call abort_physic(rname," Error: slab ocean should only be used in serial mode!",1)
	498	! endif
[1297]	499
[3100]	500	! if (is_master) write(,) trim(rname)//": Use slab-sea-ice ?"
	501	! ok_slab_sic=.true. ! default value
	502	! call getin_p("ok_slab_sic",ok_slab_sic)
	503	! if (is_master) write(,) trim(rname)//": ok_slab_sic = ",ok_slab_sic
[1297]	504
[3100]	505	! if (is_master) write(,) trim(rname)//&
	506	! ": Use heat transport for the ocean ?"
	507	! ok_slab_heat_transp=.true. ! default value
	508	! call getin_p("ok_slab_heat_transp",ok_slab_heat_transp)
	509	! if (is_master) write(,) trim(rname)//&
	510	! ": ok_slab_heat_transp = ",ok_slab_heat_transp
[1297]	511
[3299]	512	! Volcanoes (as point sources of tracers)
	513	if (is_master) write(,) trim(rname)//": Erupt volcano ?"
	514	callvolcano=.false. ! default value
	515	call getin_p("callvolcano",callvolcano)
	516	if (is_master) write(,) trim(rname)//": callvolcano = ",callvolcano
	517
	518	! Sanity check: volcano requires tracers
	519	if (callvolcano.and.(.not.tracer)) then
	520	call abort_physic(rname,&
	521	" Error: volcano option requires using tracers",1)
	522	endif
	523
[1801]	524	! Photochemistry and chemistry in the thermosphere
[1297]	525
[2583]	526	if (is_master) write(,) trim(rname)//": Use photochemistry ?"
[1801]	527	photochem=.false. ! default value
	528	call getin_p("photochem",photochem)
[2583]	529	if (is_master) write(,) trim(rname)//": photochem = ",photochem
[1297]	530
[2583]	531	if (is_master) write(,) trim(rname)//": Use photolysis heat table ?"
[2542]	532	photoheat=.false. ! default value
	533	call getin_p("photoheat",photoheat)
[2583]	534	if (is_master) write(,) trim(rname)//": photoheat = ",photoheat
[2542]	535
[2583]	536	if (is_master) write(,) trim(rname)//&
	537	": Use photolysis online calculation ?"
[2542]	538	jonline=.false. ! default value
	539	call getin_p("jonline",jonline)
[2583]	540	if (is_master) write(,) trim(rname)//": jonline = ",jonline
[2542]	541
[2583]	542	if (is_master) write(,) trim(rname)//": Use deposition ?"
[2542]	543	depos=.false. ! default value
	544	call getin_p("depos",depos)
[2583]	545	if (is_master) write(,) trim(rname)//": depos = ",depos
[2542]	546
[2583]	547	if (is_master) write(,)trim(rname)//": Production of haze ?"
[1801]	548	haze=.false. ! default value
	549	call getin_p("haze",haze)
[2583]	550	if (is_master) write(,)trim(rname)//": haze = ",haze
[1801]	551
[3309]	552	if (is_master) write(,)trim(rname)//": Output photochemical surface UV flux ?"
	553	output_writediagspecUV=.false. ! default value
	554	call getin_p("output_writediagspecUV",output_writediagspecUV)
	555	if (is_master) write(,)trim(rname)//": output_writediagspecUV = ",output_writediagspecUV
	556
[2470]	557	! Global1D mean and solar zenith angle
[1801]	558
[2470]	559	if (ngrid.eq.1) then
	560	PRINT*, 'Simulate global averaged conditions ?'
	561	global1d = .false. ! default value
	562	call getin_p("global1d",global1d)
	563	write(,) "global1d = ",global1d
	564
	565	! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
	566	if (global1d.and.diurnal) then
[2583]	567	call abort_physic(rname,'if global1d is true, diurnal must be set to false',1)
[2470]	568	endif
	569
	570	if (global1d) then
	571	PRINT *,'Solar Zenith angle (deg.) ?'
	572	PRINT *,'(assumed for averaged solar flux S/4)'
	573	szangle=60.0 ! default value
	574	call getin_p("szangle",szangle)
	575	write(,) "szangle = ",szangle
	576	endif
[2583]	577	endif ! of if (ngrid.eq.1)
[2470]	578
[253]	579	! Test of incompatibility:
	580	! if kastprof used, we must be in 1D
[1524]	581	if (kastprof.and.(ngrid.gt.1)) then
[2583]	582	call abort_physic(rname,'kastprof can only be used in 1D!',1)
[1524]	583	endif
[253]	584
[2583]	585	if (is_master) write(,)trim(rname)//&
	586	": Stratospheric temperature for kastprof mode?"
[1524]	587	Tstrat=167.0
	588	call getin_p("Tstrat",Tstrat)
[2583]	589	if (is_master) write(,)trim(rname)//": Tstrat = ",Tstrat
[253]	590
[2583]	591	if (is_master) write(,)trim(rname)//": Remove lower boundary?"
[1524]	592	nosurf=.false.
	593	call getin_p("nosurf",nosurf)
[2583]	594	if (is_master) write(,)trim(rname)//": nosurf = ",nosurf
[253]	595
	596	! Tests of incompatibility:
[1524]	597	if (nosurf.and.callsoil) then
[2583]	598	if (is_master) then
	599	write(,)trim(rname)//'nosurf not compatible with soil scheme!'
	600	write(,)trim(rname)//'... got to make a choice!'
	601	endif
	602	call abort_physic(rname,"nosurf not compatible with soil scheme!",1)
[1524]	603	endif
[253]	604
[2583]	605	if (is_master) write(,)trim(rname)//": Add an internal heat flux?", &
[1524]	606	"... matters only if callsoil=F"
	607	intheat=0.
	608	call getin_p("intheat",intheat)
[2583]	609	if (is_master) write(,)trim(rname)//": intheat = ",intheat
[952]	610
[2583]	611	if (is_master) write(,)trim(rname)//&
	612	": Use Newtonian cooling for radiative transfer?"
[1524]	613	newtonian=.false.
	614	call getin_p("newtonian",newtonian)
[2583]	615	if (is_master) write(,)trim(rname)//": newtonian = ",newtonian
[253]	616
	617	! Tests of incompatibility:
[1524]	618	if (newtonian.and.corrk) then
[2583]	619	call abort_physic(rname,'newtonian not compatible with correlated-k!',1)
[1524]	620	endif
	621	if (newtonian.and.calladj) then
[2583]	622	call abort_physic(rname,'newtonian not compatible with adjustment!',1)
[1524]	623	endif
	624	if (newtonian.and.calldifv) then
[2583]	625	call abort_physic(rname,'newtonian not compatible with a boundary layer!',1)
[1524]	626	endif
[253]	627
[2583]	628	if (is_master) write(,)trim(rname)//": Test physics timescale in 1D?"
[1524]	629	testradtimes=.false.
	630	call getin_p("testradtimes",testradtimes)
[2583]	631	if (is_master) write(,)trim(rname)//": testradtimes = ",testradtimes
[253]	632
	633	! Test of incompatibility:
	634	! if testradtimes used, we must be in 1D
[1524]	635	if (testradtimes.and.(ngrid.gt.1)) then
[2583]	636	call abort_physic(rname,'testradtimes can only be used in 1D!',1)
[1524]	637	endif
[253]	638
[2583]	639	if (is_master) write(,)trim(rname)//": Default planetary temperature?"
[1524]	640	tplanet=215.0
	641	call getin_p("tplanet",tplanet)
[2583]	642	if (is_master) write(,)trim(rname)//": tplanet = ",tplanet
[135]	643
[2583]	644	if (is_master) write(,)trim(rname)//": Which star?"
[1524]	645	startype=1 ! default value = Sol
	646	call getin_p("startype",startype)
[2583]	647	if (is_master) write(,)trim(rname)//": startype = ",startype
[135]	648
[2583]	649	if (is_master) write(,)trim(rname)//": Value of stellar flux at 1 AU?"
[1524]	650	Fat1AU=1356.0 ! default value = Sol today
	651	call getin_p("Fat1AU",Fat1AU)
[2583]	652	if (is_master) write(,)trim(rname)//": Fat1AU = ",Fat1AU
[135]	653
	654
	655	! TRACERS:
	656
[2583]	657	if (is_master) write(,)trim(rname)//": Varying H2O cloud fraction?"
[1524]	658	CLFvarying=.false. ! default value
	659	call getin_p("CLFvarying",CLFvarying)
[2583]	660	if (is_master) write(,)trim(rname)//": CLFvarying = ",CLFvarying
[253]	661
[2583]	662	if (is_master) write(,)trim(rname)//&
	663	": Value of fixed H2O cloud fraction?"
[1524]	664	CLFfixval=1.0 ! default value
	665	call getin_p("CLFfixval",CLFfixval)
[2583]	666	if (is_master) write(,)trim(rname)//": CLFfixval = ",CLFfixval
[253]	667
[2583]	668	if (is_master) write(,)trim(rname)//": fixed radii for Cloud particles?"
[1524]	669	radfixed=.false. ! default value
	670	call getin_p("radfixed",radfixed)
[2583]	671	if (is_master) write(,)trim(rname)//": radfixed = ",radfixed
[728]	672
[1524]	673	if(kastprof)then
	674	radfixed=.true.
	675	endif
[728]	676
[2583]	677	if (is_master) write(,)trim(rname)//&
	678	": Number mixing ratio of CO2 ice particles:"
[1524]	679	Nmix_co2=1.e6 ! default value
	680	call getin_p("Nmix_co2",Nmix_co2)
[2583]	681	if (is_master) write(,)trim(rname)//": Nmix_co2 = ",Nmix_co2
[135]	682
[2972]	683	if (is_master) write(,)trim(rname)//&
	684	"Number of radiatively active aerosols:"
	685	naerkind=0 ! default value
	686	call getin_p("naerkind",naerkind)
	687	if (is_master) write(,)trim(rname)//": naerkind = ",naerkind
[726]	688
[2583]	689	if (is_master) write(,)trim(rname)//": Opacity of dust (if used):"
[1524]	690	dusttau=0. ! default value
	691	call getin_p("dusttau",dusttau)
[2583]	692	if (is_master) write(,)trim(rname)//": dusttau = ",dusttau
[253]	693
[2583]	694	if (is_master) write(,)trim(rname)//": Radiatively active CO2 aerosols?"
[1524]	695	aeroco2=.false. ! default value
	696	call getin_p("aeroco2",aeroco2)
[2583]	697	if (is_master) write(,)trim(rname)//": aeroco2 = ",aeroco2
[253]	698
[2583]	699	if (is_master) write(,)trim(rname)//": Fixed CO2 aerosol distribution?"
[1524]	700	aerofixco2=.false. ! default value
	701	call getin_p("aerofixco2",aerofixco2)
[2583]	702	if (is_master) write(,)trim(rname)//": aerofixco2 = ",aerofixco2
[726]	703
[2583]	704	if (is_master) write(,)trim(rname)//": Radiatively active water ice?"
[1524]	705	aeroh2o=.false. ! default value
	706	call getin_p("aeroh2o",aeroh2o)
[2583]	707	if (is_master) write(,)trim(rname)//": aeroh2o = ",aeroh2o
[726]	708
[2583]	709	if (is_master) write(,)trim(rname)//": Fixed H2O aerosol distribution?"
[1524]	710	aerofixh2o=.false. ! default value
	711	call getin_p("aerofixh2o",aerofixh2o)
[2583]	712	if (is_master) write(,)trim(rname)//": aerofixh2o = ",aerofixh2o
[726]	713
[2583]	714	if (is_master) write(,)trim(rname)//&
	715	": Radiatively active sulfuric acid aerosols?"
[1524]	716	aeroh2so4=.false. ! default value
	717	call getin_p("aeroh2so4",aeroh2so4)
[2583]	718	if (is_master) write(,)trim(rname)//": aeroh2so4 = ",aeroh2so4
[2297]	719
[2583]	720	if (is_master) write(,)trim(rname)//&
	721	": Radiatively active auroral aerosols?"
[2297]	722	aeroaurora=.false. ! default value
	723	call getin_p("aeroaurora",aeroaurora)
[2583]	724	if (is_master) write(,)trim(rname)//": aeroaurora = ",aeroaurora
[2297]	725
[2803]	726	if (is_master) write(,)trim(rname)//&
	727	": Radiatively active Generic Condensable aerosols?"
	728	aerogeneric=0 ! default value
	729	call getin_p("aerogeneric",aerogeneric)
	730	if (is_master) write(,)trim(rname)//":aerogeneric",aerogeneric
	731
	732
[1151]	733	!=================================
[2297]	734	! TWOLAY scheme and NH3 cloudare left for retrocompatibility only,
	735	! You should now use N-LAYER scheme (see below).
[1151]	736
[2583]	737	if (is_master) write(,)trim(rname)//&
	738	": Radiatively active two-layer aerosols?"
[1524]	739	aeroback2lay=.false. ! default value
	740	call getin_p("aeroback2lay",aeroback2lay)
[2583]	741	if (is_master) write(,)trim(rname)//": aeroback2lay = ",aeroback2lay
[726]	742
[2583]	743	if (aeroback2lay.and.is_master) then
[2297]	744	print*,'Warning : The TWOLAY AEROSOL scheme is deprecated and buggy...'
	745	print*,'You should use the generic n-layer scheme (see aeronlay).'
	746	endif
	747
[2583]	748	if (is_master) write(,)trim(rname)//": Radiatively active ammonia cloud?"
[1677]	749	aeronh3=.false. ! default value
	750	call getin_p("aeronh3",aeronh3)
[2583]	751	if (is_master) write(,)trim(rname)//": aeronh3 = ",aeronh3
[1677]	752
[2583]	753	if (aeronh3.and.is_master) then
[2297]	754	print*,'Warning : You are using specific NH3 cloud scheme ...'
	755	print*,'You should use the generic n-layer scheme (see aeronlay).'
	756	endif
[1677]	757
[2831]	758	if (is_master) write(,)"Radiatively active Venus clouds ?"
	759	aerovenus=.false. ! default value
	760	call getin_p("aerovenus",aerovenus)
	761	if (aerovenus) then
	762	aerovenus1=.true. ! default value
	763	aerovenus2=.true. ! default value
	764	aerovenus2p=.true. ! default value
	765	aerovenus3=.true. ! default value
	766	aerovenusUV=.true. ! default value
	767	else
	768	aerovenus1=.false. ! default value
	769	aerovenus2=.false. ! default value
	770	aerovenus2p=.false. ! default value
	771	aerovenus3=.false. ! default value
	772	aerovenusUV=.false. ! default value
	773	endif
	774	! in case the user wants to specifically set/unset sub-options
	775	call getin_p("aerovenus1",aerovenus1)
	776	if (is_master) write(,)" aerovenus1 = ",aerovenus1
	777	call getin_p("aerovenus2",aerovenus2)
	778	if (is_master) write(,)" aerovenus2 = ",aerovenus2
	779	call getin_p("aerovenus2p",aerovenus2p)
	780	if (is_master) write(,)" aerovenus2p= ",aerovenus2p
	781	call getin_p("aerovenus3",aerovenus3)
	782	if (is_master) write(,)" aerovenus3 = ",aerovenus3
	783	call getin_p("aerovenusUV",aerovenusUV)
	784	if (is_master) write(,)" aerovenusUV= ",aerovenusUV
	785	! Sanity check: if any of the aerovenus* is set to true
	786	! then aeronovenus should also be set to true
	787	if ((.not.aerovenus).and.(aerovenus1.or.aerovenus2.or.aerovenus2p.or.&
	788	aerovenus3.or.aerovenusUV)) then
	789	if(is_master) then
	790	write(,)" Error, if you set some of the aerovenus* to true"
	791	write(,)" then flag aerovenus should be set to true as well!"
	792	endif
	793	call abort_physic("inifis"," aerovenus* flags mismatch!",1)
	794	endif
	795
[2583]	796	if (is_master) write(,)trim(rname)//&
	797	": TWOLAY AEROSOL: total optical depth "//&
	798	"in the tropospheric layer (visible)"
[1524]	799	obs_tau_col_tropo=8.D0
	800	call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
[2583]	801	if (is_master) write(,)trim(rname)//&
	802	": obs_tau_col_tropo = ",obs_tau_col_tropo
[1026]	803
[2583]	804	if (is_master) write(,)trim(rname)//&
	805	": TWOLAY AEROSOL: total optical depth "//&
	806	"in the stratospheric layer (visible)"
[1524]	807	obs_tau_col_strato=0.08D0
	808	call getin_p("obs_tau_col_strato",obs_tau_col_strato)
[2583]	809	if (is_master) write(,)trim(rname)//&
	810	": obs_tau_col_strato = ",obs_tau_col_strato
[1151]	811
[2583]	812	if (is_master) write(,)trim(rname)//&
	813	": TWOLAY AEROSOL: optprop_back2lay_vis?"
[2062]	814	optprop_back2lay_vis = 'optprop_saturn_vis_n20.dat'
	815	call getin_p("optprop_back2lay_vis",optprop_back2lay_vis)
[2583]	816	if (is_master) write(,)trim(rname)//&
	817	": optprop_back2lay_vis = ",trim(optprop_back2lay_vis)
[2062]	818
[2583]	819	if (is_master) write(,)trim(rname)//&
	820	": TWOLAY AEROSOL: optprop_back2lay_ir?"
[2062]	821	optprop_back2lay_ir = 'optprop_saturn_ir_n20.dat'
	822	call getin_p("optprop_back2lay_ir",optprop_back2lay_ir)
[2583]	823	if (is_master) write(,)trim(rname)//&
	824	": optprop_back2lay_ir = ",trim(optprop_back2lay_ir)
[2062]	825
[2583]	826	if (is_master) write(,)trim(rname)//&
	827	": TWOLAY AEROSOL: pres_bottom_tropo? in pa"
[1524]	828	pres_bottom_tropo=66000.0
	829	call getin_p("pres_bottom_tropo",pres_bottom_tropo)
[2583]	830	if (is_master) write(,)trim(rname)//&
	831	": pres_bottom_tropo = ",pres_bottom_tropo
[1151]	832
[2583]	833	if (is_master) write(,)trim(rname)//&
	834	": TWOLAY AEROSOL: pres_top_tropo? in pa"
[1524]	835	pres_top_tropo=18000.0
	836	call getin_p("pres_top_tropo",pres_top_tropo)
[2583]	837	if (is_master) write(,)trim(rname)//&
	838	": pres_top_tropo = ",pres_top_tropo
[1151]	839
[2583]	840	if (is_master) write(,)trim(rname)//&
	841	": TWOLAY AEROSOL: pres_bottom_strato? in pa"
[1524]	842	pres_bottom_strato=2000.0
	843	call getin_p("pres_bottom_strato",pres_bottom_strato)
[2583]	844	if (is_master) write(,)trim(rname)//&
	845	": pres_bottom_strato = ",pres_bottom_strato
[1151]	846
[2254]	847	! Sanity check
	848	if (pres_bottom_strato .gt. pres_top_tropo) then
[2583]	849	if(is_master) then
[2254]	850	print*,'Error : TWOLAY AEROSOL, Please ensure that in callphys.def'
	851	print*,'you have pres_top_tropo > pres_bottom_strato !'
[2583]	852	endif
	853	call abort_physic(rname," pres_top_tropo > pres_bottom_strato!",1)
[2254]	854	endif
	855
[2583]	856	if (is_master) write(,)trim(rname)//&
	857	": TWOLAY AEROSOL: pres_top_strato? in pa"
[1524]	858	pres_top_strato=100.0
	859	call getin_p("pres_top_strato",pres_top_strato)
[2583]	860	if (is_master) write(,)trim(rname)//&
	861	": pres_top_strato = ",pres_top_strato
[1151]	862
[2583]	863	if (is_master) write(,)trim(rname)//&
	864	": TWOLAY AEROSOL: particle size in the ", &
	865	"tropospheric layer, in meters"
[1524]	866	size_tropo=2.e-6
	867	call getin_p("size_tropo",size_tropo)
[2583]	868	if (is_master) write(,)trim(rname)//": size_tropo = ",size_tropo
[1151]	869
[2583]	870	if (is_master) write(,)trim(rname)//&
	871	": TWOLAY AEROSOL: particle size in the ", &
	872	"stratospheric layer, in meters"
[1524]	873	size_strato=1.e-7
	874	call getin_p("size_strato",size_strato)
[2583]	875	if (is_master) write(,)trim(rname)//": size_strato = ",size_strato
[1151]	876
[2583]	877	if (is_master) write(,)trim(rname)//&
	878	": NH3 (thin) cloud: total optical depth"
[1677]	879	tau_nh3_cloud=7.D0
	880	call getin_p("tau_nh3_cloud",tau_nh3_cloud)
[2583]	881	if (is_master) write(,)trim(rname)//": tau_nh3_cloud = ",tau_nh3_cloud
[1677]	882
[2583]	883	if (is_master) write(,)trim(rname)//": NH3 (thin) cloud pressure level"
[1677]	884	pres_nh3_cloud=7.D0
	885	call getin_p("pres_nh3_cloud",pres_nh3_cloud)
[2583]	886	if (is_master) write(,)trim(rname)//": pres_nh3_cloud = ",pres_nh3_cloud
[1677]	887
[2583]	888	if (is_master) write(,)trim(rname)//": NH3 (thin) cloud: particle sizes"
[1677]	889	size_nh3_cloud=3.e-6
	890	call getin_p("size_nh3_cloud",size_nh3_cloud)
[2583]	891	if (is_master) write(,)trim(rname)//": size_nh3_cloud = ",size_nh3_cloud
[1677]	892
[1151]	893	!=================================
[2297]	894	! Generic N-LAYER aerosol scheme
[1151]	895
[2583]	896	if (is_master) write(,)trim(rname)//&
	897	": Radiatively active generic n-layer aerosols?"
[2297]	898	aeronlay=.false. ! default value
	899	call getin_p("aeronlay",aeronlay)
[2583]	900	if (is_master) write(,)trim(rname)//": aeronlay = ",aeronlay
[2297]	901
[2583]	902	if (is_master) write(,)trim(rname)//&
	903	": Number of generic aerosols layers?"
[2297]	904	nlayaero=1 ! default value
	905	call getin_p("nlayaero",nlayaero)
	906	! Avoid to allocate arrays of size 0
	907	if (aeronlay .and. nlayaero.lt.1) then
[2583]	908	if (is_master) then
[2297]	909	print*, " You are trying to set no generic aerosols..."
	910	print*, " Set aeronlay=.false. instead ! I abort."
[2583]	911	endif
	912	call abort_physic(rname,"no generic aerosols...",1)
[2297]	913	endif
	914	if (.not. aeronlay) nlayaero=1
[2583]	915	if (is_master) write(,)trim(rname)//": nlayaero = ",nlayaero
[2297]	916
	917	! This is necessary, we just set the number of aerosol layers
	918	IF(.NOT.ALLOCATED(aeronlay_tauref)) ALLOCATE(aeronlay_tauref(nlayaero))
	919	IF(.NOT.ALLOCATED(aeronlay_lamref)) ALLOCATE(aeronlay_lamref(nlayaero))
	920	IF(.NOT.ALLOCATED(aeronlay_choice)) ALLOCATE(aeronlay_choice(nlayaero))
	921	IF(.NOT.ALLOCATED(aeronlay_pbot)) ALLOCATE(aeronlay_pbot(nlayaero))
	922	IF(.NOT.ALLOCATED(aeronlay_ptop)) ALLOCATE(aeronlay_ptop(nlayaero))
	923	IF(.NOT.ALLOCATED(aeronlay_sclhght)) ALLOCATE(aeronlay_sclhght(nlayaero))
	924	IF(.NOT.ALLOCATED(aeronlay_size)) ALLOCATE(aeronlay_size(nlayaero))
[2340]	925	IF(.NOT.ALLOCATED(aeronlay_nueff)) ALLOCATE(aeronlay_nueff(nlayaero))
[2297]	926	IF(.NOT.ALLOCATED(optprop_aeronlay_ir)) ALLOCATE(optprop_aeronlay_ir(nlayaero))
	927	IF(.NOT.ALLOCATED(optprop_aeronlay_vis)) ALLOCATE(optprop_aeronlay_vis(nlayaero))
	928
[2583]	929	if (is_master) write(,)trim(rname)//&
	930	": Generic n-layer aerosols: Optical depth at reference wavelenght"
[2297]	931	aeronlay_tauref=1.0E-1
	932	call getin_p("aeronlay_tauref",aeronlay_tauref)
[2583]	933	if (is_master) write(,)trim(rname)//&
	934	": aeronlay_tauref = ",aeronlay_tauref
[2297]	935
[2583]	936	if (is_master) write(,)trim(rname)//&
	937	": Generic n-layer aerosols: Reference wavelenght for optical depths (m)"
[2297]	938	aeronlay_lamref=0.6E-6
	939	call getin_p("aeronlay_lamref",aeronlay_lamref)
[2583]	940	if (is_master) write(,)trim(rname)//&
	941	": aeronlay_lamref = ",aeronlay_lamref
[2297]	942
[2583]	943	if (is_master) then
	944	write(,)trim(rname)//&
	945	": Generic n-layer aerosols: Vertical profile choice : "
	946	write(,)trim(rname)//&
	947	" (1) Tau btwn ptop and pbot follows atm. scale height"
	948	write(,)trim(rname)//&
	949	" or (2) Tau above pbot follows its own scale height"
	950	endif
[2297]	951	aeronlay_choice=1
	952	call getin_p("aeronlay_choice",aeronlay_choice)
[2583]	953	if (is_master) write(,)trim(rname)//&
	954	": aeronlay_choice = ",aeronlay_choice
[2297]	955
[2583]	956	if (is_master) write(,)trim(rname)//&
	957	": Generic n-layer aerosols: bottom pressures (Pa)"
[2297]	958	aeronlay_pbot=2000.0
	959	call getin_p("aeronlay_pbot",aeronlay_pbot)
[2583]	960	if (is_master) write(,)trim(rname)//": aeronlay_pbot = ",aeronlay_pbot
[2297]	961
[2583]	962	if (is_master) write(,)trim(rname)//&
	963	": Generic n-layer aerosols: (if choice=1) Top pressures (Pa) "
[2297]	964	aeronlay_ptop=300000.0
	965	call getin_p("aeronlay_ptop",aeronlay_ptop)
[2583]	966	if (is_master) write(,)trim(rname)//": aeronlay_ptop = ",aeronlay_ptop
[2297]	967
[2583]	968	if (is_master) write(,)trim(rname)//&
	969	": Generic n-layer aerosols: (if choice=2) Scale height / atm. scale height"
[2309]	970	aeronlay_sclhght=0.2
[2297]	971	call getin_p("aeronlay_sclhght",aeronlay_sclhght)
[2583]	972	if (is_master) write(,)trim(rname)//&
	973	": aeronlay_sclhght = ",aeronlay_sclhght
[2297]	974
[2583]	975	if (is_master) write(,)trim(rname)//&
	976	": Generic n-layer aerosols: particles effective radii (m)"
[2297]	977	aeronlay_size=1.e-6
	978	call getin_p("aeronlay_size",aeronlay_size)
[2583]	979	if (is_master) write(,)trim(rname)//": aeronlay_size = ",aeronlay_size
[2297]	980
[2583]	981	if (is_master) write(,)trim(rname)//&
	982	": Generic n-layer aerosols: particles radii effective variance"
[2340]	983	aeronlay_nueff=0.1
	984	call getin_p("aeronlay_nueff",aeronlay_nueff)
[2583]	985	if (is_master) write(,)trim(rname)//": aeronlay_nueff = ",aeronlay_nueff
[2340]	986
[2583]	987	if (is_master) write(,)trim(rname)//&
	988	": Generic n-layer aerosols: VIS optical properties file"
[2297]	989	optprop_aeronlay_vis = 'optprop_saturn_vis_n20.dat'
	990	call getin_p("optprop_aeronlay_vis",optprop_aeronlay_vis)
[2583]	991	if (is_master) write(,)trim(rname)//&
	992	": optprop_aeronlay_vis = ",optprop_aeronlay_vis
[2297]	993
[2583]	994	if (is_master) write(,)trim(rname)//&
	995	": Generic n-layer aerosols: IR optical properties file"
[2297]	996	optprop_aeronlay_ir = 'optprop_saturn_ir_n20.dat'
	997	call getin_p("optprop_aeronlay_ir",optprop_aeronlay_ir)
[2583]	998	if (is_master) write(,)trim(rname)//&
	999	": optprop_aeronlay_ir = ",optprop_aeronlay_ir
[2297]	1000
	1001
	1002	!=================================
	1003
[2583]	1004	if (is_master) write(,)trim(rname)//&
	1005	": Cloud pressure level (with kastprof only):"
[1524]	1006	cloudlvl=0. ! default value
	1007	call getin_p("cloudlvl",cloudlvl)
[2583]	1008	if (is_master) write(,)trim(rname)//": cloudlvl = ",cloudlvl
[305]	1009
[2583]	1010	if (is_master) write(,)trim(rname)//&
	1011	": Is the variable gas species radiatively active?"
[1524]	1012	Tstrat=167.0
	1013	varactive=.false.
	1014	call getin_p("varactive",varactive)
[2583]	1015	if (is_master) write(,)trim(rname)//": varactive = ",varactive
[135]	1016
[2583]	1017	if (is_master) write(,)trim(rname)//&
	1018	": Is the variable gas species distribution set?"
[1524]	1019	varfixed=.false.
	1020	call getin_p("varfixed",varfixed)
[2583]	1021	if (is_master) write(,)trim(rname)//": varfixed = ",varfixed
[135]	1022
[2583]	1023	if (is_master) write(,)trim(rname)//&
	1024	": What is the saturation % of the variable species?"
[1524]	1025	satval=0.8
	1026	call getin_p("satval",satval)
[2583]	1027	if (is_master) write(,)trim(rname)//": satval = ",satval
[135]	1028
	1029
	1030	! Test of incompatibility:
	1031	! if varactive, then varfixed should be false
[1524]	1032	if (varactive.and.varfixed) then
[2583]	1033	call abort_physic(rname,'if varactive, varfixed must be OFF!',1)
[1524]	1034	endif
[135]	1035
[2583]	1036	if (is_master) write(,)trim(rname)//": Gravitationnal sedimentation ?"
[1524]	1037	sedimentation=.false. ! default value
	1038	call getin_p("sedimentation",sedimentation)
[2583]	1039	if (is_master) write(,)trim(rname)//": sedimentation = ",sedimentation
[135]	1040
[2698]	1041	if (is_master) write(,)trim(rname)//": Generic Condensation of tracers ?"
	1042	generic_condensation=.false. !default value
	1043	call getin_p("generic_condensation",generic_condensation)
	1044	if (is_master) write(,)trim(rname)//": generic_condensation = ",generic_condensation
[2721]	1045
[3277]	1046	if (is_master) write(,)trim(rname)//": Generic rain of tracers ?"
[2721]	1047	generic_rain=.false. !default value
	1048	call getin_p("generic_rain",generic_rain)
	1049	if (is_master) write(,)trim(rname)//": generic_rain = ",generic_rain
[3277]	1050
	1051	if (is_master) write(,)trim(rname)//": Include moist adjustement for GCS ?"
	1052	moistadjustment_generic=.false. ! default value
	1053	call getin_p("moistadjustment_generic",moistadjustment_generic)
	1054	if (is_master) write(,)trim(rname)//": moistadjustment_generic = ", moistadjustment_generic
[3279]	1055
	1056	if (is_master) write(,)trim(rname)//": Moist convection inhibition for GCS ?"
	1057	moist_convection_inhibition=.false. ! default value
	1058	call getin_p("moist_convection_inhibition",moist_convection_inhibition)
	1059	if (is_master) write(,)trim(rname)//": moist_convection_inhibition = ", moist_convection_inhibition
[3236]	1060
	1061	if (is_master) write(,)trim(rname)//": Virtual correction ?"
	1062	virtual_correction=.false. !default value
	1063	call getin_p("virtual_correction",virtual_correction)
	1064	if (is_master) write(,)trim(rname)//": virtual_correction = ",virtual_correction
[2721]	1065
[2583]	1066	if (is_master) write(,)trim(rname)//": Compute water cycle ?"
[1524]	1067	water=.false. ! default value
	1068	call getin_p("water",water)
[2583]	1069	if (is_master) write(,)trim(rname)//": water = ",water
[135]	1070
[622]	1071	! Test of incompatibility:
	1072	! if water is true, there should be at least a tracer
[1524]	1073	if (water.and.(.not.tracer)) then
[2583]	1074	call abort_physic(rname,'if water is ON, tracer must be ON too!',1)
[1524]	1075	endif
[622]	1076
[2583]	1077	if (is_master) write(,)trim(rname)//": Include water condensation ?"
[1524]	1078	watercond=.false. ! default value
	1079	call getin_p("watercond",watercond)
[2583]	1080	if (is_master) write(,)trim(rname)//": watercond = ",watercond
[135]	1081
[622]	1082	! Test of incompatibility:
	1083	! if watercond is used, then water should be used too
[1524]	1084	if (watercond.and.(.not.water)) then
[2583]	1085	call abort_physic(rname,&
	1086	'if watercond is used, water should be used too',1)
[1524]	1087	endif
[622]	1088
[3087]	1089	if (is_master) write(,)trim(rname)//": Include water moist adjustement ?"
	1090	moistadjustment=.true. ! default value
[3099]	1091	call getin_p("moistadjustment",moistadjustment)
[3087]	1092	if (is_master) write(,)trim(rname)//": moistadjustment = ", moistadjustment
	1093
[2583]	1094	if (is_master) write(,)trim(rname)//": Include water precipitation ?"
[1524]	1095	waterrain=.false. ! default value
	1096	call getin_p("waterrain",waterrain)
[2583]	1097	if (is_master) write(,)trim(rname)//": waterrain = ",waterrain
[135]	1098
[2583]	1099	if (is_master) write(,)trim(rname)//": Include surface hydrology ?"
[1524]	1100	hydrology=.false. ! default value
	1101	call getin_p("hydrology",hydrology)
[2583]	1102	if (is_master) write(,)trim(rname)//": hydrology = ",hydrology
[135]	1103
[2583]	1104	if (is_master) write(,)trim(rname)//": Spectral Dependant albedo ?"
[1524]	1105	albedo_spectral_mode=.false. ! default value
	1106	call getin_p("albedo_spectral_mode",albedo_spectral_mode)
[2583]	1107	if (is_master) write(,)trim(rname)//&
	1108	": albedo_spectral_mode = ",albedo_spectral_mode
[486]	1109
[2583]	1110	if (is_master) then
	1111	write(,)trim(rname)//": Snow albedo ?"
	1112	write(,)trim(rname)//": If albedo_spectral_mode=.true., then this "
	1113	write(,)trim(rname)//": corresponds to the 0.5 microns snow albedo."
	1114	endif
[1524]	1115	albedosnow=0.5 ! default value
	1116	call getin_p("albedosnow",albedosnow)
[2583]	1117	if (is_master) write(,)trim(rname)//": albedosnow = ",albedosnow
[1482]	1118
[2583]	1119	if (is_master) write(,)trim(rname)//": Ocean albedo ?"
[2482]	1120	alb_ocean=0.07 ! default value
	1121	call getin_p("alb_ocean",alb_ocean)
[2583]	1122	if (is_master) write(,)trim(rname)//": alb_ocean = ",alb_ocean
[2482]	1123
[2583]	1124	if (is_master) write(,)trim(rname)//": CO2 ice albedo ?"
[1524]	1125	albedoco2ice=0.5 ! default value
	1126	call getin_p("albedoco2ice",albedoco2ice)
[2583]	1127	if (is_master) write(,)trim(rname)//": albedoco2ice = ",albedoco2ice
[135]	1128
[2583]	1129	if (is_master) write(,)trim(rname)//": Maximum ice thickness ?"
[1524]	1130	maxicethick=2.0 ! default value
	1131	call getin_p("maxicethick",maxicethick)
[2583]	1132	if (is_master) write(,)trim(rname)//": maxicethick = ",maxicethick
[135]	1133
[2583]	1134	if (is_master) write(,)trim(rname)//": Freezing point of seawater ?"
[1524]	1135	Tsaldiff=-1.8 ! default value
	1136	call getin_p("Tsaldiff",Tsaldiff)
[2583]	1137	if (is_master) write(,)trim(rname)//": Tsaldiff = ",Tsaldiff
[135]	1138
[2583]	1139	if (is_master) write(,)trim(rname)//": Minimum eddy mix coeff in 1D ?"
[2482]	1140	kmixmin=1.0e-2 ! default value
	1141	call getin_p("kmixmin",kmixmin)
[2583]	1142	if (is_master) write(,)trim(rname)//": kmixmin = ",kmixmin
[3233]	1143
	1144	varspec=.false. !default value
	1145	call getin_p("varspec",varspec)
	1146	if (varspec) then
	1147	call getin_p("varspec_data",varspec_data)
	1148	call getin_p("nvarlayer",nvarlayer)
	1149	endif
[2482]	1150
[2635]	1151	if (is_master) write(,)'Predefined Cp from dynamics is ',cpp,'J kg^-1 K^-1'
	1152	if (is_master) write(,)'Predefined Mg from dynamics is ',mugaz,'amu'
	1153
[1524]	1154	force_cpp=.false. ! default value
	1155	call getin_p("force_cpp",force_cpp)
[2635]	1156	if (force_cpp) then
	1157	if (is_master) write(,)trim(rname)//": force_cpp = ",force_cpp
	1158	if (is_master) write(,)trim(rname)//": force_cpp is deprecated.",&
	1159	"Set cpp_mugaz_mode=1 in callfis to emulate force_cpp=.true."
	1160	call abort_physic(rname,"Anyway, you need to set force_cpp=.false. to continue.",1)
	1161	endif
[589]	1162
[2635]	1163	if (is_master) write(,)trim(rname)//&
	1164	": where do you want your cpp/mugaz value to come from?",&
	1165	"=> 0: dynamics (3d), 1: forced in callfis (1d), 2: computed from gases.def (1d)?"
	1166	cpp_mugaz_mode = 0 ! default value
	1167	call getin_p("cpp_mugaz_mode",cpp_mugaz_mode)
	1168	if (is_master) write(,)trim(rname)//": cpp_mugaz_mode = ",cpp_mugaz_mode
	1169
	1170	if ((cpp_mugaz_mode >= 1).and.(is_master).and.(ngrid>1)) then
	1171	write(,)' !!!!! Be aware that having different values for cpp and mugaz in the dynamics and physics'
	1172	write(,)' in 3D can result in very pathological behavior. You have been warned !!!!!'
[2664]	1173	else if ((cpp_mugaz_mode == 0).and.(is_master).and.(ngrid==1)) then
	1174	! for this specific 1D error we will remove run.def before aborting JL22
	1175	call system("rm -rf run.def")
	1176	call abort_physic(rname,"cpp_mugaz_mode must be >= 1 in 1d",1)
[2635]	1177	endif
	1178
	1179	if (cpp_mugaz_mode == 1) then
[1524]	1180	mugaz = -99999.
[2583]	1181	if (is_master) write(,)trim(rname)//&
	1182	": MEAN MOLECULAR MASS in g mol-1 ?"
[1524]	1183	call getin_p("mugaz",mugaz)
	1184	IF (mugaz.eq.-99999.) THEN
[2635]	1185	call abort_physic(rname,"mugaz must be set if cpp_mugaz_mode = 1",1)
[1524]	1186	ENDIF
	1187	cpp = -99999.
[2583]	1188	if (is_master) write(,)trim(rname)//&
	1189	": SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
[1524]	1190	call getin_p("cpp",cpp)
	1191	IF (cpp.eq.-99999.) THEN
[2635]	1192	PRINT *, "cpp must be set if cpp_mugaz_mode = 1"
[1524]	1193	STOP
	1194	ENDIF
[2635]	1195	if (is_master) write(,)'New Cp from callfis is ',cpp,'J kg^-1 K^-1'
	1196	if (is_master) write(,)'New Mg from callfis is ',mugaz,'amu'
	1197
	1198	endif ! of if (cpp_mugaz_mode == 1)
[1524]	1199	call su_gases
	1200	call calc_cpp_mugaz
[135]	1201
[2583]	1202	if (is_master) then
	1203	PRINT*,'--------------------------------------------'
	1204	PRINT*
	1205	PRINT*
	1206	endif
[1524]	1207	ELSE
[2583]	1208	call abort_physic(rname,'Cannot read file callphys.def. Is it here ?',1)
[1524]	1209	ENDIF ! of IF(iscallphys)
[135]	1210
[2583]	1211	if (is_master) then
	1212	PRINT*
	1213	PRINT*,'inifis: daysec',daysec
	1214	PRINT*
	1215	PRINT*,'inifis: The radiative transfer is computed:'
	1216	PRINT*,' each ',iradia,' physical time-step'
	1217	PRINT,' or each ',iradiadtphys,' seconds'
	1218	PRINT*
	1219	endif
[135]	1220
	1221	!-----------------------------------------------------------------------
	1222	! Some more initialization:
	1223	! ------------------------
	1224
[1542]	1225	! Initializations for comgeomfi_h
[1832]	1226	#ifndef MESOSCALE
[1542]	1227	totarea=SSUM(ngrid,parea,1)
	1228	call planetwide_sumval(parea,totarea_planet)
[787]	1229
[1524]	1230	!! those are defined in comdiurn_h.F90
	1231	IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
	1232	IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
	1233	IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
	1234	IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
[787]	1235
[1524]	1236	DO ig=1,ngrid
	1237	sinlat(ig)=sin(plat(ig))
	1238	coslat(ig)=cos(plat(ig))
	1239	sinlon(ig)=sin(plon(ig))
	1240	coslon(ig)=cos(plon(ig))
	1241	ENDDO
[1832]	1242	#endif
[1529]	1243	! initialize variables in radinc_h
[2283]	1244	call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
[2972]	1245
	1246	! initialize variables and allocate arrays in radcommon_h
	1247	call ini_radcommon_h(naerkind)
	1248
[1524]	1249	! allocate "comsoil_h" arrays
	1250	call ini_comsoil_h(ngrid)
[1801]	1251
[1524]	1252	END SUBROUTINE inifis
[135]	1253
[1524]	1254	END MODULE inifis_mod

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

Download in other formats:

Original Format