Context Navigation

inifis_mod.F90 @ 2613

Last change on this file since 2613 was 2595, checked in by emillour, 3 years ago

Generic GCM:
Fixes and improvements in the Non-orographic GW scheme, namely:

increments are not tendencies
missing rho factor in EP flux computation
missing rho at launch altitude
changed inputs, because R and Cp are needed to compute rho and BV

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

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

Download in other formats:

Original Format