Context Navigation

inifis_mod.F90 @ 3267

Last change on this file since 3267 was 3236, checked in by gmilcareck, 2 years ago
Add virtual correction for convective adjustment and turbulent diffusion (turbdiff and vdifc) + correction of allocated tables in physiq_mod for varspec.
File size: 48.0 KB

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

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

Download in other formats:

Original Format