Context Navigation

inifis_mod.F90 @ 3073

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

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

Download in other formats:

Original Format