Context Navigation

inifis_mod.F90 @ 2803

Last change on this file since 2803 was 2803, checked in by aslmd, 2 years ago

Initialisation of Radiative Generic Condensable Aerosols

We can activate the scheme by putting aerogeneric = # of aerosols in callphys.def. This is the only needed thing for activating the radiative effects. They must be tracer in modern tracer.def

Commented out the abort if we use more than 4 aerosols

Added reading of optprop files for Radiative Generic Condensable Aerosols

We use the same file for IR and VI channel. For now, only MnS, Cr,Fe and Mg2SiO4 can be read. If you want to add another specie, check the code, it is explained how to quickly do that (right above the Initializations)

Added radii calculation for Radiative Generic Condensable aerosols

Changed the hardcoded size of the totalemit array

The hardcoded size is now 1900 instead of 100 so we don't exceed the array size when working at high spectral resolution (very rare case)

Added opacity computation for Radiative Generic Condensable aerosols

We do this computation in the same fashion as what's performed on water and dust.

switch iniaerosol and initracer order, to prepare for the RGCS scheme

Needed to switch the order of initialization so we can use the RGCS scheme without the assumption that ice and vap tracer of the same specie are following each other in the traceur.def file

File size: 44.2 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)//": call radiative transfer ?"
[1524]	243	callrad=.true. ! default value
	244	call getin_p("callrad",callrad)
[2583]	245	if (is_master) write(,) trim(rname)//": callrad = ",callrad
[135]	246
[2583]	247	if (is_master) write(,) trim(rname)//&
	248	": call correlated-k radiative transfer ?"
[1524]	249	corrk=.true. ! default value
	250	call getin_p("corrk",corrk)
[2583]	251	if (is_master) write(,) trim(rname)//": corrk = ",corrk
[135]	252
[2583]	253	if (is_master) write(,) trim(rname)//&
	254	": prohibit calculations outside corrk T grid?"
[1524]	255	strictboundcorrk=.true. ! default value
	256	call getin_p("strictboundcorrk",strictboundcorrk)
[2583]	257	if (is_master) write(,) trim(rname)//&
	258	": strictboundcorrk = ",strictboundcorrk
[1145]	259
[2583]	260	if (is_master) write(,) trim(rname)//&
[2655]	261	": prohibit calculations outside CIA T grid?"
	262	strictboundcia=.true. ! default value
	263	call getin_p("strictboundcia",strictboundcia)
[2633]	264	if (is_master) write(,) trim(rname)//&
[2655]	265	": strictboundcia = ",strictboundcia
[2633]	266
	267	if (is_master) write(,) trim(rname)//&
[2583]	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)
[2655]	306	endif
	307
	308	if (is_master) write(,) trim(rname)//&
	309	": CIA - normal or equilibrium ortho-para mixture for H2?"
	310	H2orthopara_mixture = 'normal'
	311	call getin_p("H2orthopara_mixture",H2orthopara_mixture)
	312	if (is_master) write(,)trim(rname)//&
	313	": H2orthopara_mixture = ",trim(H2orthopara_mixture)
[2245]	314
[2583]	315	if (is_master) write(,) trim(rname)//&
	316	": call turbulent vertical diffusion ?"
[1524]	317	calldifv=.true. ! default value
	318	call getin_p("calldifv",calldifv)
[2583]	319	if (is_master) write(,) trim(rname)//": calldifv = ",calldifv
[135]	320
[2583]	321	if (is_master) write(,) trim(rname)//": use turbdiff instead of vdifc ?"
[1524]	322	UseTurbDiff=.true. ! default value
	323	call getin_p("UseTurbDiff",UseTurbDiff)
[2583]	324	if (is_master) write(,) trim(rname)//": UseTurbDiff = ",UseTurbDiff
[596]	325
[2583]	326	if (is_master) write(,) trim(rname)//": call convective adjustment ?"
[1524]	327	calladj=.true. ! default value
	328	call getin_p("calladj",calladj)
[2583]	329	if (is_master) write(,) trim(rname)//": calladj = ",calladj
[135]	330
[2583]	331	if (is_master) write(,)trim(rname)//&
	332	": call Lott's non-oro GWs parameterisation scheme ?"
[2299]	333	calllott_nonoro=.false. ! default value
	334	call getin_p("calllott_nonoro",calllott_nonoro)
[2583]	335	if (is_master) write(,)trim(rname)//&
	336	": calllott_nonoro = ",calllott_nonoro
[2299]	337
[2583]	338	if (is_master) write(,) trim(rname)//": call thermal plume model ?"
[2060]	339	calltherm=.false. ! default value
	340	call getin_p("calltherm",calltherm)
[2583]	341	if (is_master) write(,) trim(rname)//": calltherm = ",calltherm
[2060]	342
[2583]	343	if (is_master) write(,) trim(rname)//": call CO2 condensation ?"
[1524]	344	co2cond=.false. ! default value
	345	call getin_p("co2cond",co2cond)
[2583]	346	if (is_master) write(,) trim(rname)//": co2cond = ",co2cond
[650]	347	! Test of incompatibility
[1524]	348	if (co2cond.and.(.not.tracer)) then
[2583]	349	call abort_physic(rname,"Error: We need a CO2 ice tracer to condense CO2!",1)
[1524]	350	endif
[716]	351
[2583]	352	if (is_master) write(,) trim(rname)//": CO2 supersaturation level ?"
[1524]	353	co2supsat=1.0 ! default value
	354	call getin_p("co2supsat",co2supsat)
[2583]	355	if (is_master) write(,) trim(rname)//": co2supsat = ",co2supsat
[253]	356
[2583]	357	if (is_master) write(,) trim(rname)//&
	358	": Radiative timescale for Newtonian cooling (in Earth days)?"
[1524]	359	tau_relax=30. ! default value
	360	call getin_p("tau_relax",tau_relax)
[2583]	361	if (is_master) write(,) trim(rname)//": tau_relax = ",tau_relax
[1524]	362	tau_relax=tau_relax243600 ! convert Earth days --> seconds
[253]	363
[2583]	364	if (is_master) write(,)trim(rname)//&
	365	": call thermal conduction in the soil ?"
[1524]	366	callsoil=.true. ! default value
	367	call getin_p("callsoil",callsoil)
[2583]	368	if (is_master) write(,) trim(rname)//": callsoil = ",callsoil
[135]	369
[2583]	370	if (is_master) write(,)trim(rname)//&
	371	": Rad transfer is computed every iradia", &
	372	" physical timestep"
[1524]	373	iradia=1 ! default value
	374	call getin_p("iradia",iradia)
[2583]	375	if (is_master) write(,)trim(rname)//": iradia = ",iradia
[135]	376
[2583]	377	if (is_master) write(,)trim(rname)//": Rayleigh scattering ?"
[1524]	378	rayleigh=.false.
	379	call getin_p("rayleigh",rayleigh)
[2583]	380	if (is_master) write(,)trim(rname)//": rayleigh = ",rayleigh
[135]	381
[2583]	382	if (is_master) write(,) trim(rname)//&
	383	": Use blackbody for stellar spectrum ?"
[1524]	384	stelbbody=.false. ! default value
	385	call getin_p("stelbbody",stelbbody)
[2583]	386	if (is_master) write(,) trim(rname)//": stelbbody = ",stelbbody
[135]	387
[2583]	388	if (is_master) write(,) trim(rname)//": Stellar blackbody temperature ?"
[1524]	389	stelTbb=5800.0 ! default value
	390	call getin_p("stelTbb",stelTbb)
[2583]	391	if (is_master) write(,) trim(rname)//": stelTbb = ",stelTbb
[253]	392
[2583]	393	if (is_master) write(,) trim(rname)//": Output mean OLR in 1D?"
[1524]	394	meanOLR=.false.
	395	call getin_p("meanOLR",meanOLR)
[2583]	396	if (is_master) write(,) trim(rname)//": meanOLR = ",meanOLR
[135]	397
[2583]	398	if (is_master) write(,)trim(rname)//": Output spectral OLR in 3D?"
[1524]	399	specOLR=.false.
	400	call getin_p("specOLR",specOLR)
[2583]	401	if (is_master) write(,)trim(rname)//": specOLR = ",specOLR
[135]	402
[2583]	403	if (is_master) write(,)trim(rname)//&
	404	": Output diagnostic optical thickness attenuation exp(-tau)?"
[2131]	405	diagdtau=.false.
	406	call getin_p("diagdtau",diagdtau)
[2583]	407	if (is_master) write(,)trim(rname)//": diagdtau = ",diagdtau
[2131]	408
[2583]	409	if (is_master) write(,)trim(rname)//": Operate in kastprof mode?"
[1524]	410	kastprof=.false.
	411	call getin_p("kastprof",kastprof)
[2583]	412	if (is_master) write(,)trim(rname)//": kastprof = ",kastprof
[253]	413
[2583]	414	if (is_master) write(,)trim(rname)//&
	415	": Uniform absorption in radiative transfer?"
[1524]	416	graybody=.false.
	417	call getin_p("graybody",graybody)
[2583]	418	if (is_master) write(,)trim(rname)//": graybody = ",graybody
[253]	419
[1538]	420	! Soil model
[2583]	421	if (is_master) write(,)trim(rname)//&
	422	": Number of sub-surface layers for soil scheme?"
[1538]	423	! default value of nsoilmx set in comsoil_h
	424	call getin_p("nsoilmx",nsoilmx)
[2583]	425	if (is_master) write(,)trim(rname)//": nsoilmx=",nsoilmx
[1538]	426
[2583]	427	if (is_master) write(,)trim(rname)//&
	428	": Thickness of topmost soil layer (m)?"
[1538]	429	! default value of lay1_soil set in comsoil_h
	430	call getin_p("lay1_soil",lay1_soil)
[2583]	431	if (is_master) write(,)trim(rname)//": lay1_soil = ",lay1_soil
[1538]	432
[2583]	433	if (is_master) write(,)trim(rname)//&
	434	": Coefficient for soil layer thickness distribution?"
[1538]	435	! default value of alpha_soil set in comsoil_h
	436	call getin_p("alpha_soil",alpha_soil)
[2583]	437	if (is_master) write(,)trim(rname)//": alpha_soil = ",alpha_soil
[1538]	438
[1297]	439	! Slab Ocean
[2583]	440	if (is_master) write(,) trim(rname)//": Use slab-ocean ?"
[1524]	441	ok_slab_ocean=.false. ! default value
	442	call getin_p("ok_slab_ocean",ok_slab_ocean)
[2583]	443	if (is_master) write(,) trim(rname)//": ok_slab_ocean = ",ok_slab_ocean
[1529]	444	! Sanity check: for now slab oncean only works in serial mode
	445	if (ok_slab_ocean.and.is_parallel) then
[2583]	446	call abort_physic(rname," Error: slab ocean should only be used in serial mode!",1)
[1529]	447	endif
[1297]	448
[2583]	449	if (is_master) write(,) trim(rname)//": Use slab-sea-ice ?"
[1524]	450	ok_slab_sic=.true. ! default value
	451	call getin_p("ok_slab_sic",ok_slab_sic)
[2583]	452	if (is_master) write(,) trim(rname)//": ok_slab_sic = ",ok_slab_sic
[1297]	453
[2583]	454	if (is_master) write(,) trim(rname)//&
	455	": Use heat transport for the ocean ?"
[1524]	456	ok_slab_heat_transp=.true. ! default value
	457	call getin_p("ok_slab_heat_transp",ok_slab_heat_transp)
[2583]	458	if (is_master) write(,) trim(rname)//&
	459	": ok_slab_heat_transp = ",ok_slab_heat_transp
[1297]	460
[1801]	461	! Photochemistry and chemistry in the thermosphere
[1297]	462
[2583]	463	if (is_master) write(,) trim(rname)//": Use photochemistry ?"
[1801]	464	photochem=.false. ! default value
	465	call getin_p("photochem",photochem)
[2583]	466	if (is_master) write(,) trim(rname)//": photochem = ",photochem
[1297]	467
[2583]	468	if (is_master) write(,) trim(rname)//": Use photolysis heat table ?"
[2542]	469	photoheat=.false. ! default value
	470	call getin_p("photoheat",photoheat)
[2583]	471	if (is_master) write(,) trim(rname)//": photoheat = ",photoheat
[2542]	472
[2583]	473	if (is_master) write(,) trim(rname)//&
	474	": Use photolysis online calculation ?"
[2542]	475	jonline=.false. ! default value
	476	call getin_p("jonline",jonline)
[2583]	477	if (is_master) write(,) trim(rname)//": jonline = ",jonline
[2542]	478
[2583]	479	if (is_master) write(,) trim(rname)//": Use deposition ?"
[2542]	480	depos=.false. ! default value
	481	call getin_p("depos",depos)
[2583]	482	if (is_master) write(,) trim(rname)//": depos = ",depos
[2542]	483
[2583]	484	if (is_master) write(,)trim(rname)//": Production of haze ?"
[1801]	485	haze=.false. ! default value
	486	call getin_p("haze",haze)
[2583]	487	if (is_master) write(,)trim(rname)//": haze = ",haze
[1801]	488
[2470]	489	! Global1D mean and solar zenith angle
[1801]	490
[2470]	491	if (ngrid.eq.1) then
	492	PRINT*, 'Simulate global averaged conditions ?'
	493	global1d = .false. ! default value
	494	call getin_p("global1d",global1d)
	495	write(,) "global1d = ",global1d
	496
	497	! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
	498	if (global1d.and.diurnal) then
[2583]	499	call abort_physic(rname,'if global1d is true, diurnal must be set to false',1)
[2470]	500	endif
	501
	502	if (global1d) then
	503	PRINT *,'Solar Zenith angle (deg.) ?'
	504	PRINT *,'(assumed for averaged solar flux S/4)'
	505	szangle=60.0 ! default value
	506	call getin_p("szangle",szangle)
	507	write(,) "szangle = ",szangle
	508	endif
[2583]	509	endif ! of if (ngrid.eq.1)
[2470]	510
[253]	511	! Test of incompatibility:
	512	! if kastprof used, we must be in 1D
[1524]	513	if (kastprof.and.(ngrid.gt.1)) then
[2583]	514	call abort_physic(rname,'kastprof can only be used in 1D!',1)
[1524]	515	endif
[253]	516
[2583]	517	if (is_master) write(,)trim(rname)//&
	518	": Stratospheric temperature for kastprof mode?"
[1524]	519	Tstrat=167.0
	520	call getin_p("Tstrat",Tstrat)
[2583]	521	if (is_master) write(,)trim(rname)//": Tstrat = ",Tstrat
[253]	522
[2583]	523	if (is_master) write(,)trim(rname)//": Remove lower boundary?"
[1524]	524	nosurf=.false.
	525	call getin_p("nosurf",nosurf)
[2583]	526	if (is_master) write(,)trim(rname)//": nosurf = ",nosurf
[253]	527
	528	! Tests of incompatibility:
[1524]	529	if (nosurf.and.callsoil) then
[2583]	530	if (is_master) then
	531	write(,)trim(rname)//'nosurf not compatible with soil scheme!'
	532	write(,)trim(rname)//'... got to make a choice!'
	533	endif
	534	call abort_physic(rname,"nosurf not compatible with soil scheme!",1)
[1524]	535	endif
[253]	536
[2583]	537	if (is_master) write(,)trim(rname)//": Add an internal heat flux?", &
[1524]	538	"... matters only if callsoil=F"
	539	intheat=0.
	540	call getin_p("intheat",intheat)
[2583]	541	if (is_master) write(,)trim(rname)//": intheat = ",intheat
[952]	542
[2583]	543	if (is_master) write(,)trim(rname)//&
	544	": Use Newtonian cooling for radiative transfer?"
[1524]	545	newtonian=.false.
	546	call getin_p("newtonian",newtonian)
[2583]	547	if (is_master) write(,)trim(rname)//": newtonian = ",newtonian
[253]	548
	549	! Tests of incompatibility:
[1524]	550	if (newtonian.and.corrk) then
[2583]	551	call abort_physic(rname,'newtonian not compatible with correlated-k!',1)
[1524]	552	endif
	553	if (newtonian.and.calladj) then
[2583]	554	call abort_physic(rname,'newtonian not compatible with adjustment!',1)
[1524]	555	endif
	556	if (newtonian.and.calldifv) then
[2583]	557	call abort_physic(rname,'newtonian not compatible with a boundary layer!',1)
[1524]	558	endif
[253]	559
[2583]	560	if (is_master) write(,)trim(rname)//": Test physics timescale in 1D?"
[1524]	561	testradtimes=.false.
	562	call getin_p("testradtimes",testradtimes)
[2583]	563	if (is_master) write(,)trim(rname)//": testradtimes = ",testradtimes
[253]	564
	565	! Test of incompatibility:
	566	! if testradtimes used, we must be in 1D
[1524]	567	if (testradtimes.and.(ngrid.gt.1)) then
[2583]	568	call abort_physic(rname,'testradtimes can only be used in 1D!',1)
[1524]	569	endif
[253]	570
[2583]	571	if (is_master) write(,)trim(rname)//": Default planetary temperature?"
[1524]	572	tplanet=215.0
	573	call getin_p("tplanet",tplanet)
[2583]	574	if (is_master) write(,)trim(rname)//": tplanet = ",tplanet
[135]	575
[2583]	576	if (is_master) write(,)trim(rname)//": Which star?"
[1524]	577	startype=1 ! default value = Sol
	578	call getin_p("startype",startype)
[2583]	579	if (is_master) write(,)trim(rname)//": startype = ",startype
[135]	580
[2583]	581	if (is_master) write(,)trim(rname)//": Value of stellar flux at 1 AU?"
[1524]	582	Fat1AU=1356.0 ! default value = Sol today
	583	call getin_p("Fat1AU",Fat1AU)
[2583]	584	if (is_master) write(,)trim(rname)//": Fat1AU = ",Fat1AU
[135]	585
	586
	587	! TRACERS:
	588
[2583]	589	if (is_master) write(,)trim(rname)//": Varying H2O cloud fraction?"
[1524]	590	CLFvarying=.false. ! default value
	591	call getin_p("CLFvarying",CLFvarying)
[2583]	592	if (is_master) write(,)trim(rname)//": CLFvarying = ",CLFvarying
[253]	593
[2583]	594	if (is_master) write(,)trim(rname)//&
	595	": Value of fixed H2O cloud fraction?"
[1524]	596	CLFfixval=1.0 ! default value
	597	call getin_p("CLFfixval",CLFfixval)
[2583]	598	if (is_master) write(,)trim(rname)//": CLFfixval = ",CLFfixval
[253]	599
[2583]	600	if (is_master) write(,)trim(rname)//": fixed radii for Cloud particles?"
[1524]	601	radfixed=.false. ! default value
	602	call getin_p("radfixed",radfixed)
[2583]	603	if (is_master) write(,)trim(rname)//": radfixed = ",radfixed
[728]	604
[1524]	605	if(kastprof)then
	606	radfixed=.true.
	607	endif
[728]	608
[2583]	609	if (is_master) write(,)trim(rname)//&
	610	": Number mixing ratio of CO2 ice particles:"
[1524]	611	Nmix_co2=1.e6 ! default value
	612	call getin_p("Nmix_co2",Nmix_co2)
[2583]	613	if (is_master) write(,)trim(rname)//": Nmix_co2 = ",Nmix_co2
[135]	614
[726]	615	! write(,)"Number of radiatively active aerosols:"
	616	! naerkind=0. ! default value
[1315]	617	! call getin_p("naerkind",naerkind)
[726]	618	! write(,)" naerkind = ",naerkind
	619
[2583]	620	if (is_master) write(,)trim(rname)//": Opacity of dust (if used):"
[1524]	621	dusttau=0. ! default value
	622	call getin_p("dusttau",dusttau)
[2583]	623	if (is_master) write(,)trim(rname)//": dusttau = ",dusttau
[253]	624
[2583]	625	if (is_master) write(,)trim(rname)//": Radiatively active CO2 aerosols?"
[1524]	626	aeroco2=.false. ! default value
	627	call getin_p("aeroco2",aeroco2)
[2583]	628	if (is_master) write(,)trim(rname)//": aeroco2 = ",aeroco2
[253]	629
[2583]	630	if (is_master) write(,)trim(rname)//": Fixed CO2 aerosol distribution?"
[1524]	631	aerofixco2=.false. ! default value
	632	call getin_p("aerofixco2",aerofixco2)
[2583]	633	if (is_master) write(,)trim(rname)//": aerofixco2 = ",aerofixco2
[726]	634
[2583]	635	if (is_master) write(,)trim(rname)//": Radiatively active water ice?"
[1524]	636	aeroh2o=.false. ! default value
	637	call getin_p("aeroh2o",aeroh2o)
[2583]	638	if (is_master) write(,)trim(rname)//": aeroh2o = ",aeroh2o
[726]	639
[2583]	640	if (is_master) write(,)trim(rname)//": Fixed H2O aerosol distribution?"
[1524]	641	aerofixh2o=.false. ! default value
	642	call getin_p("aerofixh2o",aerofixh2o)
[2583]	643	if (is_master) write(,)trim(rname)//": aerofixh2o = ",aerofixh2o
[726]	644
[2583]	645	if (is_master) write(,)trim(rname)//&
	646	": Radiatively active sulfuric acid aerosols?"
[1524]	647	aeroh2so4=.false. ! default value
	648	call getin_p("aeroh2so4",aeroh2so4)
[2583]	649	if (is_master) write(,)trim(rname)//": aeroh2so4 = ",aeroh2so4
[2297]	650
[2583]	651	if (is_master) write(,)trim(rname)//&
	652	": Radiatively active auroral aerosols?"
[2297]	653	aeroaurora=.false. ! default value
	654	call getin_p("aeroaurora",aeroaurora)
[2583]	655	if (is_master) write(,)trim(rname)//": aeroaurora = ",aeroaurora
[2297]	656
[2803]	657	if (is_master) write(,)trim(rname)//&
	658	": Radiatively active Generic Condensable aerosols?"
	659	aerogeneric=0 ! default value
	660	call getin_p("aerogeneric",aerogeneric)
	661	if (is_master) write(,)trim(rname)//":aerogeneric",aerogeneric
	662
	663
[1151]	664	!=================================
[2297]	665	! TWOLAY scheme and NH3 cloudare left for retrocompatibility only,
	666	! You should now use N-LAYER scheme (see below).
[1151]	667
[2583]	668	if (is_master) write(,)trim(rname)//&
	669	": Radiatively active two-layer aerosols?"
[1524]	670	aeroback2lay=.false. ! default value
	671	call getin_p("aeroback2lay",aeroback2lay)
[2583]	672	if (is_master) write(,)trim(rname)//": aeroback2lay = ",aeroback2lay
[726]	673
[2583]	674	if (aeroback2lay.and.is_master) then
[2297]	675	print*,'Warning : The TWOLAY AEROSOL scheme is deprecated and buggy...'
	676	print*,'You should use the generic n-layer scheme (see aeronlay).'
	677	endif
	678
[2583]	679	if (is_master) write(,)trim(rname)//": Radiatively active ammonia cloud?"
[1677]	680	aeronh3=.false. ! default value
	681	call getin_p("aeronh3",aeronh3)
[2583]	682	if (is_master) write(,)trim(rname)//": aeronh3 = ",aeronh3
[1677]	683
[2583]	684	if (aeronh3.and.is_master) then
[2297]	685	print*,'Warning : You are using specific NH3 cloud scheme ...'
	686	print*,'You should use the generic n-layer scheme (see aeronlay).'
	687	endif
[1677]	688
[2583]	689	if (is_master) write(,)trim(rname)//&
	690	": TWOLAY AEROSOL: total optical depth "//&
	691	"in the tropospheric layer (visible)"
[1524]	692	obs_tau_col_tropo=8.D0
	693	call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
[2583]	694	if (is_master) write(,)trim(rname)//&
	695	": obs_tau_col_tropo = ",obs_tau_col_tropo
[1026]	696
[2583]	697	if (is_master) write(,)trim(rname)//&
	698	": TWOLAY AEROSOL: total optical depth "//&
	699	"in the stratospheric layer (visible)"
[1524]	700	obs_tau_col_strato=0.08D0
	701	call getin_p("obs_tau_col_strato",obs_tau_col_strato)
[2583]	702	if (is_master) write(,)trim(rname)//&
	703	": obs_tau_col_strato = ",obs_tau_col_strato
[1151]	704
[2583]	705	if (is_master) write(,)trim(rname)//&
	706	": TWOLAY AEROSOL: optprop_back2lay_vis?"
[2062]	707	optprop_back2lay_vis = 'optprop_saturn_vis_n20.dat'
	708	call getin_p("optprop_back2lay_vis",optprop_back2lay_vis)
[2583]	709	if (is_master) write(,)trim(rname)//&
	710	": optprop_back2lay_vis = ",trim(optprop_back2lay_vis)
[2062]	711
[2583]	712	if (is_master) write(,)trim(rname)//&
	713	": TWOLAY AEROSOL: optprop_back2lay_ir?"
[2062]	714	optprop_back2lay_ir = 'optprop_saturn_ir_n20.dat'
	715	call getin_p("optprop_back2lay_ir",optprop_back2lay_ir)
[2583]	716	if (is_master) write(,)trim(rname)//&
	717	": optprop_back2lay_ir = ",trim(optprop_back2lay_ir)
[2062]	718
[2583]	719	if (is_master) write(,)trim(rname)//&
	720	": TWOLAY AEROSOL: pres_bottom_tropo? in pa"
[1524]	721	pres_bottom_tropo=66000.0
	722	call getin_p("pres_bottom_tropo",pres_bottom_tropo)
[2583]	723	if (is_master) write(,)trim(rname)//&
	724	": pres_bottom_tropo = ",pres_bottom_tropo
[1151]	725
[2583]	726	if (is_master) write(,)trim(rname)//&
	727	": TWOLAY AEROSOL: pres_top_tropo? in pa"
[1524]	728	pres_top_tropo=18000.0
	729	call getin_p("pres_top_tropo",pres_top_tropo)
[2583]	730	if (is_master) write(,)trim(rname)//&
	731	": pres_top_tropo = ",pres_top_tropo
[1151]	732
[2583]	733	if (is_master) write(,)trim(rname)//&
	734	": TWOLAY AEROSOL: pres_bottom_strato? in pa"
[1524]	735	pres_bottom_strato=2000.0
	736	call getin_p("pres_bottom_strato",pres_bottom_strato)
[2583]	737	if (is_master) write(,)trim(rname)//&
	738	": pres_bottom_strato = ",pres_bottom_strato
[1151]	739
[2254]	740	! Sanity check
	741	if (pres_bottom_strato .gt. pres_top_tropo) then
[2583]	742	if(is_master) then
[2254]	743	print*,'Error : TWOLAY AEROSOL, Please ensure that in callphys.def'
	744	print*,'you have pres_top_tropo > pres_bottom_strato !'
[2583]	745	endif
	746	call abort_physic(rname," pres_top_tropo > pres_bottom_strato!",1)
[2254]	747	endif
	748
[2583]	749	if (is_master) write(,)trim(rname)//&
	750	": TWOLAY AEROSOL: pres_top_strato? in pa"
[1524]	751	pres_top_strato=100.0
	752	call getin_p("pres_top_strato",pres_top_strato)
[2583]	753	if (is_master) write(,)trim(rname)//&
	754	": pres_top_strato = ",pres_top_strato
[1151]	755
[2583]	756	if (is_master) write(,)trim(rname)//&
	757	": TWOLAY AEROSOL: particle size in the ", &
	758	"tropospheric layer, in meters"
[1524]	759	size_tropo=2.e-6
	760	call getin_p("size_tropo",size_tropo)
[2583]	761	if (is_master) write(,)trim(rname)//": size_tropo = ",size_tropo
[1151]	762
[2583]	763	if (is_master) write(,)trim(rname)//&
	764	": TWOLAY AEROSOL: particle size in the ", &
	765	"stratospheric layer, in meters"
[1524]	766	size_strato=1.e-7
	767	call getin_p("size_strato",size_strato)
[2583]	768	if (is_master) write(,)trim(rname)//": size_strato = ",size_strato
[1151]	769
[2583]	770	if (is_master) write(,)trim(rname)//&
	771	": NH3 (thin) cloud: total optical depth"
[1677]	772	tau_nh3_cloud=7.D0
	773	call getin_p("tau_nh3_cloud",tau_nh3_cloud)
[2583]	774	if (is_master) write(,)trim(rname)//": tau_nh3_cloud = ",tau_nh3_cloud
[1677]	775
[2583]	776	if (is_master) write(,)trim(rname)//": NH3 (thin) cloud pressure level"
[1677]	777	pres_nh3_cloud=7.D0
	778	call getin_p("pres_nh3_cloud",pres_nh3_cloud)
[2583]	779	if (is_master) write(,)trim(rname)//": pres_nh3_cloud = ",pres_nh3_cloud
[1677]	780
[2583]	781	if (is_master) write(,)trim(rname)//": NH3 (thin) cloud: particle sizes"
[1677]	782	size_nh3_cloud=3.e-6
	783	call getin_p("size_nh3_cloud",size_nh3_cloud)
[2583]	784	if (is_master) write(,)trim(rname)//": size_nh3_cloud = ",size_nh3_cloud
[1677]	785
[1151]	786	!=================================
[2297]	787	! Generic N-LAYER aerosol scheme
[1151]	788
[2583]	789	if (is_master) write(,)trim(rname)//&
	790	": Radiatively active generic n-layer aerosols?"
[2297]	791	aeronlay=.false. ! default value
	792	call getin_p("aeronlay",aeronlay)
[2583]	793	if (is_master) write(,)trim(rname)//": aeronlay = ",aeronlay
[2297]	794
[2583]	795	if (is_master) write(,)trim(rname)//&
	796	": Number of generic aerosols layers?"
[2297]	797	nlayaero=1 ! default value
	798	call getin_p("nlayaero",nlayaero)
	799	! Avoid to allocate arrays of size 0
	800	if (aeronlay .and. nlayaero.lt.1) then
[2583]	801	if (is_master) then
[2297]	802	print*, " You are trying to set no generic aerosols..."
	803	print*, " Set aeronlay=.false. instead ! I abort."
[2583]	804	endif
	805	call abort_physic(rname,"no generic aerosols...",1)
[2297]	806	endif
	807	if (.not. aeronlay) nlayaero=1
[2583]	808	if (is_master) write(,)trim(rname)//": nlayaero = ",nlayaero
[2297]	809
	810	! This is necessary, we just set the number of aerosol layers
	811	IF(.NOT.ALLOCATED(aeronlay_tauref)) ALLOCATE(aeronlay_tauref(nlayaero))
	812	IF(.NOT.ALLOCATED(aeronlay_lamref)) ALLOCATE(aeronlay_lamref(nlayaero))
	813	IF(.NOT.ALLOCATED(aeronlay_choice)) ALLOCATE(aeronlay_choice(nlayaero))
	814	IF(.NOT.ALLOCATED(aeronlay_pbot)) ALLOCATE(aeronlay_pbot(nlayaero))
	815	IF(.NOT.ALLOCATED(aeronlay_ptop)) ALLOCATE(aeronlay_ptop(nlayaero))
	816	IF(.NOT.ALLOCATED(aeronlay_sclhght)) ALLOCATE(aeronlay_sclhght(nlayaero))
	817	IF(.NOT.ALLOCATED(aeronlay_size)) ALLOCATE(aeronlay_size(nlayaero))
[2340]	818	IF(.NOT.ALLOCATED(aeronlay_nueff)) ALLOCATE(aeronlay_nueff(nlayaero))
[2297]	819	IF(.NOT.ALLOCATED(optprop_aeronlay_ir)) ALLOCATE(optprop_aeronlay_ir(nlayaero))
	820	IF(.NOT.ALLOCATED(optprop_aeronlay_vis)) ALLOCATE(optprop_aeronlay_vis(nlayaero))
	821
[2583]	822	if (is_master) write(,)trim(rname)//&
	823	": Generic n-layer aerosols: Optical depth at reference wavelenght"
[2297]	824	aeronlay_tauref=1.0E-1
	825	call getin_p("aeronlay_tauref",aeronlay_tauref)
[2583]	826	if (is_master) write(,)trim(rname)//&
	827	": aeronlay_tauref = ",aeronlay_tauref
[2297]	828
[2583]	829	if (is_master) write(,)trim(rname)//&
	830	": Generic n-layer aerosols: Reference wavelenght for optical depths (m)"
[2297]	831	aeronlay_lamref=0.6E-6
	832	call getin_p("aeronlay_lamref",aeronlay_lamref)
[2583]	833	if (is_master) write(,)trim(rname)//&
	834	": aeronlay_lamref = ",aeronlay_lamref
[2297]	835
[2583]	836	if (is_master) then
	837	write(,)trim(rname)//&
	838	": Generic n-layer aerosols: Vertical profile choice : "
	839	write(,)trim(rname)//&
	840	" (1) Tau btwn ptop and pbot follows atm. scale height"
	841	write(,)trim(rname)//&
	842	" or (2) Tau above pbot follows its own scale height"
	843	endif
[2297]	844	aeronlay_choice=1
	845	call getin_p("aeronlay_choice",aeronlay_choice)
[2583]	846	if (is_master) write(,)trim(rname)//&
	847	": aeronlay_choice = ",aeronlay_choice
[2297]	848
[2583]	849	if (is_master) write(,)trim(rname)//&
	850	": Generic n-layer aerosols: bottom pressures (Pa)"
[2297]	851	aeronlay_pbot=2000.0
	852	call getin_p("aeronlay_pbot",aeronlay_pbot)
[2583]	853	if (is_master) write(,)trim(rname)//": aeronlay_pbot = ",aeronlay_pbot
[2297]	854
[2583]	855	if (is_master) write(,)trim(rname)//&
	856	": Generic n-layer aerosols: (if choice=1) Top pressures (Pa) "
[2297]	857	aeronlay_ptop=300000.0
	858	call getin_p("aeronlay_ptop",aeronlay_ptop)
[2583]	859	if (is_master) write(,)trim(rname)//": aeronlay_ptop = ",aeronlay_ptop
[2297]	860
[2583]	861	if (is_master) write(,)trim(rname)//&
	862	": Generic n-layer aerosols: (if choice=2) Scale height / atm. scale height"
[2309]	863	aeronlay_sclhght=0.2
[2297]	864	call getin_p("aeronlay_sclhght",aeronlay_sclhght)
[2583]	865	if (is_master) write(,)trim(rname)//&
	866	": aeronlay_sclhght = ",aeronlay_sclhght
[2297]	867
[2583]	868	if (is_master) write(,)trim(rname)//&
	869	": Generic n-layer aerosols: particles effective radii (m)"
[2297]	870	aeronlay_size=1.e-6
	871	call getin_p("aeronlay_size",aeronlay_size)
[2583]	872	if (is_master) write(,)trim(rname)//": aeronlay_size = ",aeronlay_size
[2297]	873
[2583]	874	if (is_master) write(,)trim(rname)//&
	875	": Generic n-layer aerosols: particles radii effective variance"
[2340]	876	aeronlay_nueff=0.1
	877	call getin_p("aeronlay_nueff",aeronlay_nueff)
[2583]	878	if (is_master) write(,)trim(rname)//": aeronlay_nueff = ",aeronlay_nueff
[2340]	879
[2583]	880	if (is_master) write(,)trim(rname)//&
	881	": Generic n-layer aerosols: VIS optical properties file"
[2297]	882	optprop_aeronlay_vis = 'optprop_saturn_vis_n20.dat'
	883	call getin_p("optprop_aeronlay_vis",optprop_aeronlay_vis)
[2583]	884	if (is_master) write(,)trim(rname)//&
	885	": optprop_aeronlay_vis = ",optprop_aeronlay_vis
[2297]	886
[2583]	887	if (is_master) write(,)trim(rname)//&
	888	": Generic n-layer aerosols: IR optical properties file"
[2297]	889	optprop_aeronlay_ir = 'optprop_saturn_ir_n20.dat'
	890	call getin_p("optprop_aeronlay_ir",optprop_aeronlay_ir)
[2583]	891	if (is_master) write(,)trim(rname)//&
	892	": optprop_aeronlay_ir = ",optprop_aeronlay_ir
[2297]	893
	894
	895	!=================================
	896
[2583]	897	if (is_master) write(,)trim(rname)//&
	898	": Cloud pressure level (with kastprof only):"
[1524]	899	cloudlvl=0. ! default value
	900	call getin_p("cloudlvl",cloudlvl)
[2583]	901	if (is_master) write(,)trim(rname)//": cloudlvl = ",cloudlvl
[305]	902
[2583]	903	if (is_master) write(,)trim(rname)//&
	904	": Is the variable gas species radiatively active?"
[1524]	905	Tstrat=167.0
	906	varactive=.false.
	907	call getin_p("varactive",varactive)
[2583]	908	if (is_master) write(,)trim(rname)//": varactive = ",varactive
[135]	909
[2583]	910	if (is_master) write(,)trim(rname)//&
	911	": Is the variable gas species distribution set?"
[1524]	912	varfixed=.false.
	913	call getin_p("varfixed",varfixed)
[2583]	914	if (is_master) write(,)trim(rname)//": varfixed = ",varfixed
[135]	915
[2583]	916	if (is_master) write(,)trim(rname)//&
	917	": What is the saturation % of the variable species?"
[1524]	918	satval=0.8
	919	call getin_p("satval",satval)
[2583]	920	if (is_master) write(,)trim(rname)//": satval = ",satval
[135]	921
	922
	923	! Test of incompatibility:
	924	! if varactive, then varfixed should be false
[1524]	925	if (varactive.and.varfixed) then
[2583]	926	call abort_physic(rname,'if varactive, varfixed must be OFF!',1)
[1524]	927	endif
[135]	928
[2583]	929	if (is_master) write(,)trim(rname)//": Gravitationnal sedimentation ?"
[1524]	930	sedimentation=.false. ! default value
	931	call getin_p("sedimentation",sedimentation)
[2583]	932	if (is_master) write(,)trim(rname)//": sedimentation = ",sedimentation
[135]	933
[2698]	934	if (is_master) write(,)trim(rname)//": Generic Condensation of tracers ?"
	935	generic_condensation=.false. !default value
	936	call getin_p("generic_condensation",generic_condensation)
	937	if (is_master) write(,)trim(rname)//": generic_condensation = ",generic_condensation
[2721]	938
	939	if (is_master) write(,)trim(rname)//": Generic Condensation of tracers ?"
	940	generic_rain=.false. !default value
	941	call getin_p("generic_rain",generic_rain)
	942	if (is_master) write(,)trim(rname)//": generic_rain = ",generic_rain
	943
[2583]	944	if (is_master) write(,)trim(rname)//": Compute water cycle ?"
[1524]	945	water=.false. ! default value
	946	call getin_p("water",water)
[2583]	947	if (is_master) write(,)trim(rname)//": water = ",water
[135]	948
[622]	949	! Test of incompatibility:
	950	! if water is true, there should be at least a tracer
[1524]	951	if (water.and.(.not.tracer)) then
[2583]	952	call abort_physic(rname,'if water is ON, tracer must be ON too!',1)
[1524]	953	endif
[622]	954
[2583]	955	if (is_master) write(,)trim(rname)//": Include water condensation ?"
[1524]	956	watercond=.false. ! default value
	957	call getin_p("watercond",watercond)
[2583]	958	if (is_master) write(,)trim(rname)//": watercond = ",watercond
[135]	959
[622]	960	! Test of incompatibility:
	961	! if watercond is used, then water should be used too
[1524]	962	if (watercond.and.(.not.water)) then
[2583]	963	call abort_physic(rname,&
	964	'if watercond is used, water should be used too',1)
[1524]	965	endif
[622]	966
[2583]	967	if (is_master) write(,)trim(rname)//": Include water precipitation ?"
[1524]	968	waterrain=.false. ! default value
	969	call getin_p("waterrain",waterrain)
[2583]	970	if (is_master) write(,)trim(rname)//": waterrain = ",waterrain
[135]	971
[2583]	972	if (is_master) write(,)trim(rname)//": Include surface hydrology ?"
[1524]	973	hydrology=.false. ! default value
	974	call getin_p("hydrology",hydrology)
[2583]	975	if (is_master) write(,)trim(rname)//": hydrology = ",hydrology
[135]	976
[2583]	977	if (is_master) write(,)trim(rname)//": Evolve surface water sources ?"
[1524]	978	sourceevol=.false. ! default value
	979	call getin_p("sourceevol",sourceevol)
[2583]	980	if (is_master) write(,)trim(rname)//": sourceevol = ",sourceevol
[135]	981
[2583]	982	if (is_master) write(,)trim(rname)//": Ice evolution timestep ?"
[1524]	983	icetstep=100.0 ! default value
	984	call getin_p("icetstep",icetstep)
[2583]	985	if (is_master) write(,)trim(rname)//": icetstep = ",icetstep
[1498]	986
[2583]	987	if (is_master) write(,)trim(rname)//": Spectral Dependant albedo ?"
[1524]	988	albedo_spectral_mode=.false. ! default value
	989	call getin_p("albedo_spectral_mode",albedo_spectral_mode)
[2583]	990	if (is_master) write(,)trim(rname)//&
	991	": albedo_spectral_mode = ",albedo_spectral_mode
[486]	992
[2583]	993	if (is_master) then
	994	write(,)trim(rname)//": Snow albedo ?"
	995	write(,)trim(rname)//": If albedo_spectral_mode=.true., then this "
	996	write(,)trim(rname)//": corresponds to the 0.5 microns snow albedo."
	997	endif
[1524]	998	albedosnow=0.5 ! default value
	999	call getin_p("albedosnow",albedosnow)
[2583]	1000	if (is_master) write(,)trim(rname)//": albedosnow = ",albedosnow
[1482]	1001
[2583]	1002	if (is_master) write(,)trim(rname)//": Ocean albedo ?"
[2482]	1003	alb_ocean=0.07 ! default value
	1004	call getin_p("alb_ocean",alb_ocean)
[2583]	1005	if (is_master) write(,)trim(rname)//": alb_ocean = ",alb_ocean
[2482]	1006
[2583]	1007	if (is_master) write(,)trim(rname)//": CO2 ice albedo ?"
[1524]	1008	albedoco2ice=0.5 ! default value
	1009	call getin_p("albedoco2ice",albedoco2ice)
[2583]	1010	if (is_master) write(,)trim(rname)//": albedoco2ice = ",albedoco2ice
[135]	1011
[2583]	1012	if (is_master) write(,)trim(rname)//": Maximum ice thickness ?"
[1524]	1013	maxicethick=2.0 ! default value
	1014	call getin_p("maxicethick",maxicethick)
[2583]	1015	if (is_master) write(,)trim(rname)//": maxicethick = ",maxicethick
[135]	1016
[2583]	1017	if (is_master) write(,)trim(rname)//": Freezing point of seawater ?"
[1524]	1018	Tsaldiff=-1.8 ! default value
	1019	call getin_p("Tsaldiff",Tsaldiff)
[2583]	1020	if (is_master) write(,)trim(rname)//": Tsaldiff = ",Tsaldiff
[135]	1021
[2583]	1022	if (is_master) write(,)trim(rname)//": Minimum eddy mix coeff in 1D ?"
[2482]	1023	kmixmin=1.0e-2 ! default value
	1024	call getin_p("kmixmin",kmixmin)
[2583]	1025	if (is_master) write(,)trim(rname)//": kmixmin = ",kmixmin
[2482]	1026
[2635]	1027	if (is_master) write(,)'Predefined Cp from dynamics is ',cpp,'J kg^-1 K^-1'
	1028	if (is_master) write(,)'Predefined Mg from dynamics is ',mugaz,'amu'
	1029
[1524]	1030	force_cpp=.false. ! default value
	1031	call getin_p("force_cpp",force_cpp)
[2635]	1032	if (force_cpp) then
	1033	if (is_master) write(,)trim(rname)//": force_cpp = ",force_cpp
	1034	if (is_master) write(,)trim(rname)//": force_cpp is deprecated.",&
	1035	"Set cpp_mugaz_mode=1 in callfis to emulate force_cpp=.true."
	1036	call abort_physic(rname,"Anyway, you need to set force_cpp=.false. to continue.",1)
	1037	endif
[589]	1038
[2635]	1039	if (is_master) write(,)trim(rname)//&
	1040	": where do you want your cpp/mugaz value to come from?",&
	1041	"=> 0: dynamics (3d), 1: forced in callfis (1d), 2: computed from gases.def (1d)?"
	1042	cpp_mugaz_mode = 0 ! default value
	1043	call getin_p("cpp_mugaz_mode",cpp_mugaz_mode)
	1044	if (is_master) write(,)trim(rname)//": cpp_mugaz_mode = ",cpp_mugaz_mode
	1045
	1046	if ((cpp_mugaz_mode >= 1).and.(is_master).and.(ngrid>1)) then
	1047	write(,)' !!!!! Be aware that having different values for cpp and mugaz in the dynamics and physics'
	1048	write(,)' in 3D can result in very pathological behavior. You have been warned !!!!!'
[2664]	1049	else if ((cpp_mugaz_mode == 0).and.(is_master).and.(ngrid==1)) then
	1050	! for this specific 1D error we will remove run.def before aborting JL22
	1051	call system("rm -rf run.def")
	1052	call abort_physic(rname,"cpp_mugaz_mode must be >= 1 in 1d",1)
[2635]	1053	endif
	1054
	1055	if (cpp_mugaz_mode == 1) then
[1524]	1056	mugaz = -99999.
[2583]	1057	if (is_master) write(,)trim(rname)//&
	1058	": MEAN MOLECULAR MASS in g mol-1 ?"
[1524]	1059	call getin_p("mugaz",mugaz)
	1060	IF (mugaz.eq.-99999.) THEN
[2635]	1061	call abort_physic(rname,"mugaz must be set if cpp_mugaz_mode = 1",1)
[1524]	1062	ENDIF
	1063	cpp = -99999.
[2583]	1064	if (is_master) write(,)trim(rname)//&
	1065	": SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
[1524]	1066	call getin_p("cpp",cpp)
	1067	IF (cpp.eq.-99999.) THEN
[2635]	1068	PRINT *, "cpp must be set if cpp_mugaz_mode = 1"
[1524]	1069	STOP
	1070	ENDIF
[2635]	1071	if (is_master) write(,)'New Cp from callfis is ',cpp,'J kg^-1 K^-1'
	1072	if (is_master) write(,)'New Mg from callfis is ',mugaz,'amu'
	1073
	1074	endif ! of if (cpp_mugaz_mode == 1)
[1524]	1075	call su_gases
	1076	call calc_cpp_mugaz
[135]	1077
[2583]	1078	if (is_master) then
	1079	PRINT*,'--------------------------------------------'
	1080	PRINT*
	1081	PRINT*
	1082	endif
[1524]	1083	ELSE
[2583]	1084	call abort_physic(rname,'Cannot read file callphys.def. Is it here ?',1)
[1524]	1085	ENDIF ! of IF(iscallphys)
[135]	1086
[2583]	1087	if (is_master) then
	1088	PRINT*
	1089	PRINT*,'inifis: daysec',daysec
	1090	PRINT*
	1091	PRINT*,'inifis: The radiative transfer is computed:'
	1092	PRINT*,' each ',iradia,' physical time-step'
	1093	PRINT,' or each ',iradiadtphys,' seconds'
	1094	PRINT*
	1095	endif
[135]	1096
	1097	!-----------------------------------------------------------------------
	1098	! Some more initialization:
	1099	! ------------------------
	1100
[1542]	1101	! Initializations for comgeomfi_h
[1832]	1102	#ifndef MESOSCALE
[1542]	1103	totarea=SSUM(ngrid,parea,1)
	1104	call planetwide_sumval(parea,totarea_planet)
[787]	1105
[1524]	1106	!! those are defined in comdiurn_h.F90
	1107	IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
	1108	IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
	1109	IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
	1110	IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
[787]	1111
[1524]	1112	DO ig=1,ngrid
	1113	sinlat(ig)=sin(plat(ig))
	1114	coslat(ig)=cos(plat(ig))
	1115	sinlon(ig)=sin(plon(ig))
	1116	coslon(ig)=cos(plon(ig))
	1117	ENDDO
[1832]	1118	#endif
[1529]	1119	! initialize variables in radinc_h
[2283]	1120	call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
[1529]	1121
[1524]	1122	! allocate "comsoil_h" arrays
	1123	call ini_comsoil_h(ngrid)
[1801]	1124
[1524]	1125	END SUBROUTINE inifis
[135]	1126
[1524]	1127	END MODULE inifis_mod

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

Download in other formats:

Original Format