Context Navigation

inifis_mod.F90 @ 3342

Last change on this file since 3342 was 3309, checked in by yjaziri, 22 months ago

GENERIC PCM:

Cosmetic + clarifying some variables and comments in chemistry
add controle variable for actinique UV flux in photochemistry and output surface UV flux in diagspecUV.nc

File size: 49.4 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	! Volcanoes (as point sources of tracers)
506	if (is_master) write(,) trim(rname)//": Erupt volcano ?"
507	callvolcano=.false. ! default value
508	call getin_p("callvolcano",callvolcano)
509	if (is_master) write(,) trim(rname)//": callvolcano = ",callvolcano
510
511	! Sanity check: volcano requires tracers
512	if (callvolcano.and.(.not.tracer)) then
513	call abort_physic(rname,&
514	" Error: volcano option requires using tracers",1)
515	endif
516
517	! Photochemistry and chemistry in the thermosphere
518
519	if (is_master) write(,) trim(rname)//": Use photochemistry ?"
520	photochem=.false. ! default value
521	call getin_p("photochem",photochem)
522	if (is_master) write(,) trim(rname)//": photochem = ",photochem
523
524	if (is_master) write(,) trim(rname)//": Use photolysis heat table ?"
525	photoheat=.false. ! default value
526	call getin_p("photoheat",photoheat)
527	if (is_master) write(,) trim(rname)//": photoheat = ",photoheat
528
529	if (is_master) write(,) trim(rname)//&
530	": Use photolysis online calculation ?"
531	jonline=.false. ! default value
532	call getin_p("jonline",jonline)
533	if (is_master) write(,) trim(rname)//": jonline = ",jonline
534
535	if (is_master) write(,) trim(rname)//": Use deposition ?"
536	depos=.false. ! default value
537	call getin_p("depos",depos)
538	if (is_master) write(,) trim(rname)//": depos = ",depos
539
540	if (is_master) write(,)trim(rname)//": Production of haze ?"
541	haze=.false. ! default value
542	call getin_p("haze",haze)
543	if (is_master) write(,)trim(rname)//": haze = ",haze
544
545	if (is_master) write(,)trim(rname)//": Output photochemical surface UV flux ?"
546	output_writediagspecUV=.false. ! default value
547	call getin_p("output_writediagspecUV",output_writediagspecUV)
548	if (is_master) write(,)trim(rname)//": output_writediagspecUV = ",output_writediagspecUV
549
550	! Global1D mean and solar zenith angle
551
552	if (ngrid.eq.1) then
553	PRINT*, 'Simulate global averaged conditions ?'
554	global1d = .false. ! default value
555	call getin_p("global1d",global1d)
556	write(,) "global1d = ",global1d
557
558	! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
559	if (global1d.and.diurnal) then
560	call abort_physic(rname,'if global1d is true, diurnal must be set to false',1)
561	endif
562
563	if (global1d) then
564	PRINT *,'Solar Zenith angle (deg.) ?'
565	PRINT *,'(assumed for averaged solar flux S/4)'
566	szangle=60.0 ! default value
567	call getin_p("szangle",szangle)
568	write(,) "szangle = ",szangle
569	endif
570	endif ! of if (ngrid.eq.1)
571
572	! Test of incompatibility:
573	! if kastprof used, we must be in 1D
574	if (kastprof.and.(ngrid.gt.1)) then
575	call abort_physic(rname,'kastprof can only be used in 1D!',1)
576	endif
577
578	if (is_master) write(,)trim(rname)//&
579	": Stratospheric temperature for kastprof mode?"
580	Tstrat=167.0
581	call getin_p("Tstrat",Tstrat)
582	if (is_master) write(,)trim(rname)//": Tstrat = ",Tstrat
583
584	if (is_master) write(,)trim(rname)//": Remove lower boundary?"
585	nosurf=.false.
586	call getin_p("nosurf",nosurf)
587	if (is_master) write(,)trim(rname)//": nosurf = ",nosurf
588
589	! Tests of incompatibility:
590	if (nosurf.and.callsoil) then
591	if (is_master) then
592	write(,)trim(rname)//'nosurf not compatible with soil scheme!'
593	write(,)trim(rname)//'... got to make a choice!'
594	endif
595	call abort_physic(rname,"nosurf not compatible with soil scheme!",1)
596	endif
597
598	if (is_master) write(,)trim(rname)//": Add an internal heat flux?", &
599	"... matters only if callsoil=F"
600	intheat=0.
601	call getin_p("intheat",intheat)
602	if (is_master) write(,)trim(rname)//": intheat = ",intheat
603
604	if (is_master) write(,)trim(rname)//&
605	": Use Newtonian cooling for radiative transfer?"
606	newtonian=.false.
607	call getin_p("newtonian",newtonian)
608	if (is_master) write(,)trim(rname)//": newtonian = ",newtonian
609
610	! Tests of incompatibility:
611	if (newtonian.and.corrk) then
612	call abort_physic(rname,'newtonian not compatible with correlated-k!',1)
613	endif
614	if (newtonian.and.calladj) then
615	call abort_physic(rname,'newtonian not compatible with adjustment!',1)
616	endif
617	if (newtonian.and.calldifv) then
618	call abort_physic(rname,'newtonian not compatible with a boundary layer!',1)
619	endif
620
621	if (is_master) write(,)trim(rname)//": Test physics timescale in 1D?"
622	testradtimes=.false.
623	call getin_p("testradtimes",testradtimes)
624	if (is_master) write(,)trim(rname)//": testradtimes = ",testradtimes
625
626	! Test of incompatibility:
627	! if testradtimes used, we must be in 1D
628	if (testradtimes.and.(ngrid.gt.1)) then
629	call abort_physic(rname,'testradtimes can only be used in 1D!',1)
630	endif
631
632	if (is_master) write(,)trim(rname)//": Default planetary temperature?"
633	tplanet=215.0
634	call getin_p("tplanet",tplanet)
635	if (is_master) write(,)trim(rname)//": tplanet = ",tplanet
636
637	if (is_master) write(,)trim(rname)//": Which star?"
638	startype=1 ! default value = Sol
639	call getin_p("startype",startype)
640	if (is_master) write(,)trim(rname)//": startype = ",startype
641
642	if (is_master) write(,)trim(rname)//": Value of stellar flux at 1 AU?"
643	Fat1AU=1356.0 ! default value = Sol today
644	call getin_p("Fat1AU",Fat1AU)
645	if (is_master) write(,)trim(rname)//": Fat1AU = ",Fat1AU
646
647
648	! TRACERS:
649
650	if (is_master) write(,)trim(rname)//": Varying H2O cloud fraction?"
651	CLFvarying=.false. ! default value
652	call getin_p("CLFvarying",CLFvarying)
653	if (is_master) write(,)trim(rname)//": CLFvarying = ",CLFvarying
654
655	if (is_master) write(,)trim(rname)//&
656	": Value of fixed H2O cloud fraction?"
657	CLFfixval=1.0 ! default value
658	call getin_p("CLFfixval",CLFfixval)
659	if (is_master) write(,)trim(rname)//": CLFfixval = ",CLFfixval
660
661	if (is_master) write(,)trim(rname)//": fixed radii for Cloud particles?"
662	radfixed=.false. ! default value
663	call getin_p("radfixed",radfixed)
664	if (is_master) write(,)trim(rname)//": radfixed = ",radfixed
665
666	if(kastprof)then
667	radfixed=.true.
668	endif
669
670	if (is_master) write(,)trim(rname)//&
671	": Number mixing ratio of CO2 ice particles:"
672	Nmix_co2=1.e6 ! default value
673	call getin_p("Nmix_co2",Nmix_co2)
674	if (is_master) write(,)trim(rname)//": Nmix_co2 = ",Nmix_co2
675
676	if (is_master) write(,)trim(rname)//&
677	"Number of radiatively active aerosols:"
678	naerkind=0 ! default value
679	call getin_p("naerkind",naerkind)
680	if (is_master) write(,)trim(rname)//": naerkind = ",naerkind
681
682	if (is_master) write(,)trim(rname)//": Opacity of dust (if used):"
683	dusttau=0. ! default value
684	call getin_p("dusttau",dusttau)
685	if (is_master) write(,)trim(rname)//": dusttau = ",dusttau
686
687	if (is_master) write(,)trim(rname)//": Radiatively active CO2 aerosols?"
688	aeroco2=.false. ! default value
689	call getin_p("aeroco2",aeroco2)
690	if (is_master) write(,)trim(rname)//": aeroco2 = ",aeroco2
691
692	if (is_master) write(,)trim(rname)//": Fixed CO2 aerosol distribution?"
693	aerofixco2=.false. ! default value
694	call getin_p("aerofixco2",aerofixco2)
695	if (is_master) write(,)trim(rname)//": aerofixco2 = ",aerofixco2
696
697	if (is_master) write(,)trim(rname)//": Radiatively active water ice?"
698	aeroh2o=.false. ! default value
699	call getin_p("aeroh2o",aeroh2o)
700	if (is_master) write(,)trim(rname)//": aeroh2o = ",aeroh2o
701
702	if (is_master) write(,)trim(rname)//": Fixed H2O aerosol distribution?"
703	aerofixh2o=.false. ! default value
704	call getin_p("aerofixh2o",aerofixh2o)
705	if (is_master) write(,)trim(rname)//": aerofixh2o = ",aerofixh2o
706
707	if (is_master) write(,)trim(rname)//&
708	": Radiatively active sulfuric acid aerosols?"
709	aeroh2so4=.false. ! default value
710	call getin_p("aeroh2so4",aeroh2so4)
711	if (is_master) write(,)trim(rname)//": aeroh2so4 = ",aeroh2so4
712
713	if (is_master) write(,)trim(rname)//&
714	": Radiatively active auroral aerosols?"
715	aeroaurora=.false. ! default value
716	call getin_p("aeroaurora",aeroaurora)
717	if (is_master) write(,)trim(rname)//": aeroaurora = ",aeroaurora
718
719	if (is_master) write(,)trim(rname)//&
720	": Radiatively active Generic Condensable aerosols?"
721	aerogeneric=0 ! default value
722	call getin_p("aerogeneric",aerogeneric)
723	if (is_master) write(,)trim(rname)//":aerogeneric",aerogeneric
724
725
726	!=================================
727	! TWOLAY scheme and NH3 cloudare left for retrocompatibility only,
728	! You should now use N-LAYER scheme (see below).
729
730	if (is_master) write(,)trim(rname)//&
731	": Radiatively active two-layer aerosols?"
732	aeroback2lay=.false. ! default value
733	call getin_p("aeroback2lay",aeroback2lay)
734	if (is_master) write(,)trim(rname)//": aeroback2lay = ",aeroback2lay
735
736	if (aeroback2lay.and.is_master) then
737	print*,'Warning : The TWOLAY AEROSOL scheme is deprecated and buggy...'
738	print*,'You should use the generic n-layer scheme (see aeronlay).'
739	endif
740
741	if (is_master) write(,)trim(rname)//": Radiatively active ammonia cloud?"
742	aeronh3=.false. ! default value
743	call getin_p("aeronh3",aeronh3)
744	if (is_master) write(,)trim(rname)//": aeronh3 = ",aeronh3
745
746	if (aeronh3.and.is_master) then
747	print*,'Warning : You are using specific NH3 cloud scheme ...'
748	print*,'You should use the generic n-layer scheme (see aeronlay).'
749	endif
750
751	if (is_master) write(,)"Radiatively active Venus clouds ?"
752	aerovenus=.false. ! default value
753	call getin_p("aerovenus",aerovenus)
754	if (aerovenus) then
755	aerovenus1=.true. ! default value
756	aerovenus2=.true. ! default value
757	aerovenus2p=.true. ! default value
758	aerovenus3=.true. ! default value
759	aerovenusUV=.true. ! default value
760	else
761	aerovenus1=.false. ! default value
762	aerovenus2=.false. ! default value
763	aerovenus2p=.false. ! default value
764	aerovenus3=.false. ! default value
765	aerovenusUV=.false. ! default value
766	endif
767	! in case the user wants to specifically set/unset sub-options
768	call getin_p("aerovenus1",aerovenus1)
769	if (is_master) write(,)" aerovenus1 = ",aerovenus1
770	call getin_p("aerovenus2",aerovenus2)
771	if (is_master) write(,)" aerovenus2 = ",aerovenus2
772	call getin_p("aerovenus2p",aerovenus2p)
773	if (is_master) write(,)" aerovenus2p= ",aerovenus2p
774	call getin_p("aerovenus3",aerovenus3)
775	if (is_master) write(,)" aerovenus3 = ",aerovenus3
776	call getin_p("aerovenusUV",aerovenusUV)
777	if (is_master) write(,)" aerovenusUV= ",aerovenusUV
778	! Sanity check: if any of the aerovenus* is set to true
779	! then aeronovenus should also be set to true
780	if ((.not.aerovenus).and.(aerovenus1.or.aerovenus2.or.aerovenus2p.or.&
781	aerovenus3.or.aerovenusUV)) then
782	if(is_master) then
783	write(,)" Error, if you set some of the aerovenus* to true"
784	write(,)" then flag aerovenus should be set to true as well!"
785	endif
786	call abort_physic("inifis"," aerovenus* flags mismatch!",1)
787	endif
788
789	if (is_master) write(,)trim(rname)//&
790	": TWOLAY AEROSOL: total optical depth "//&
791	"in the tropospheric layer (visible)"
792	obs_tau_col_tropo=8.D0
793	call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
794	if (is_master) write(,)trim(rname)//&
795	": obs_tau_col_tropo = ",obs_tau_col_tropo
796
797	if (is_master) write(,)trim(rname)//&
798	": TWOLAY AEROSOL: total optical depth "//&
799	"in the stratospheric layer (visible)"
800	obs_tau_col_strato=0.08D0
801	call getin_p("obs_tau_col_strato",obs_tau_col_strato)
802	if (is_master) write(,)trim(rname)//&
803	": obs_tau_col_strato = ",obs_tau_col_strato
804
805	if (is_master) write(,)trim(rname)//&
806	": TWOLAY AEROSOL: optprop_back2lay_vis?"
807	optprop_back2lay_vis = 'optprop_saturn_vis_n20.dat'
808	call getin_p("optprop_back2lay_vis",optprop_back2lay_vis)
809	if (is_master) write(,)trim(rname)//&
810	": optprop_back2lay_vis = ",trim(optprop_back2lay_vis)
811
812	if (is_master) write(,)trim(rname)//&
813	": TWOLAY AEROSOL: optprop_back2lay_ir?"
814	optprop_back2lay_ir = 'optprop_saturn_ir_n20.dat'
815	call getin_p("optprop_back2lay_ir",optprop_back2lay_ir)
816	if (is_master) write(,)trim(rname)//&
817	": optprop_back2lay_ir = ",trim(optprop_back2lay_ir)
818
819	if (is_master) write(,)trim(rname)//&
820	": TWOLAY AEROSOL: pres_bottom_tropo? in pa"
821	pres_bottom_tropo=66000.0
822	call getin_p("pres_bottom_tropo",pres_bottom_tropo)
823	if (is_master) write(,)trim(rname)//&
824	": pres_bottom_tropo = ",pres_bottom_tropo
825
826	if (is_master) write(,)trim(rname)//&
827	": TWOLAY AEROSOL: pres_top_tropo? in pa"
828	pres_top_tropo=18000.0
829	call getin_p("pres_top_tropo",pres_top_tropo)
830	if (is_master) write(,)trim(rname)//&
831	": pres_top_tropo = ",pres_top_tropo
832
833	if (is_master) write(,)trim(rname)//&
834	": TWOLAY AEROSOL: pres_bottom_strato? in pa"
835	pres_bottom_strato=2000.0
836	call getin_p("pres_bottom_strato",pres_bottom_strato)
837	if (is_master) write(,)trim(rname)//&
838	": pres_bottom_strato = ",pres_bottom_strato
839
840	! Sanity check
841	if (pres_bottom_strato .gt. pres_top_tropo) then
842	if(is_master) then
843	print*,'Error : TWOLAY AEROSOL, Please ensure that in callphys.def'
844	print*,'you have pres_top_tropo > pres_bottom_strato !'
845	endif
846	call abort_physic(rname," pres_top_tropo > pres_bottom_strato!",1)
847	endif
848
849	if (is_master) write(,)trim(rname)//&
850	": TWOLAY AEROSOL: pres_top_strato? in pa"
851	pres_top_strato=100.0
852	call getin_p("pres_top_strato",pres_top_strato)
853	if (is_master) write(,)trim(rname)//&
854	": pres_top_strato = ",pres_top_strato
855
856	if (is_master) write(,)trim(rname)//&
857	": TWOLAY AEROSOL: particle size in the ", &
858	"tropospheric layer, in meters"
859	size_tropo=2.e-6
860	call getin_p("size_tropo",size_tropo)
861	if (is_master) write(,)trim(rname)//": size_tropo = ",size_tropo
862
863	if (is_master) write(,)trim(rname)//&
864	": TWOLAY AEROSOL: particle size in the ", &
865	"stratospheric layer, in meters"
866	size_strato=1.e-7
867	call getin_p("size_strato",size_strato)
868	if (is_master) write(,)trim(rname)//": size_strato = ",size_strato
869
870	if (is_master) write(,)trim(rname)//&
871	": NH3 (thin) cloud: total optical depth"
872	tau_nh3_cloud=7.D0
873	call getin_p("tau_nh3_cloud",tau_nh3_cloud)
874	if (is_master) write(,)trim(rname)//": tau_nh3_cloud = ",tau_nh3_cloud
875
876	if (is_master) write(,)trim(rname)//": NH3 (thin) cloud pressure level"
877	pres_nh3_cloud=7.D0
878	call getin_p("pres_nh3_cloud",pres_nh3_cloud)
879	if (is_master) write(,)trim(rname)//": pres_nh3_cloud = ",pres_nh3_cloud
880
881	if (is_master) write(,)trim(rname)//": NH3 (thin) cloud: particle sizes"
882	size_nh3_cloud=3.e-6
883	call getin_p("size_nh3_cloud",size_nh3_cloud)
884	if (is_master) write(,)trim(rname)//": size_nh3_cloud = ",size_nh3_cloud
885
886	!=================================
887	! Generic N-LAYER aerosol scheme
888
889	if (is_master) write(,)trim(rname)//&
890	": Radiatively active generic n-layer aerosols?"
891	aeronlay=.false. ! default value
892	call getin_p("aeronlay",aeronlay)
893	if (is_master) write(,)trim(rname)//": aeronlay = ",aeronlay
894
895	if (is_master) write(,)trim(rname)//&
896	": Number of generic aerosols layers?"
897	nlayaero=1 ! default value
898	call getin_p("nlayaero",nlayaero)
899	! Avoid to allocate arrays of size 0
900	if (aeronlay .and. nlayaero.lt.1) then
901	if (is_master) then
902	print*, " You are trying to set no generic aerosols..."
903	print*, " Set aeronlay=.false. instead ! I abort."
904	endif
905	call abort_physic(rname,"no generic aerosols...",1)
906	endif
907	if (.not. aeronlay) nlayaero=1
908	if (is_master) write(,)trim(rname)//": nlayaero = ",nlayaero
909
910	! This is necessary, we just set the number of aerosol layers
911	IF(.NOT.ALLOCATED(aeronlay_tauref)) ALLOCATE(aeronlay_tauref(nlayaero))
912	IF(.NOT.ALLOCATED(aeronlay_lamref)) ALLOCATE(aeronlay_lamref(nlayaero))
913	IF(.NOT.ALLOCATED(aeronlay_choice)) ALLOCATE(aeronlay_choice(nlayaero))
914	IF(.NOT.ALLOCATED(aeronlay_pbot)) ALLOCATE(aeronlay_pbot(nlayaero))
915	IF(.NOT.ALLOCATED(aeronlay_ptop)) ALLOCATE(aeronlay_ptop(nlayaero))
916	IF(.NOT.ALLOCATED(aeronlay_sclhght)) ALLOCATE(aeronlay_sclhght(nlayaero))
917	IF(.NOT.ALLOCATED(aeronlay_size)) ALLOCATE(aeronlay_size(nlayaero))
918	IF(.NOT.ALLOCATED(aeronlay_nueff)) ALLOCATE(aeronlay_nueff(nlayaero))
919	IF(.NOT.ALLOCATED(optprop_aeronlay_ir)) ALLOCATE(optprop_aeronlay_ir(nlayaero))
920	IF(.NOT.ALLOCATED(optprop_aeronlay_vis)) ALLOCATE(optprop_aeronlay_vis(nlayaero))
921
922	if (is_master) write(,)trim(rname)//&
923	": Generic n-layer aerosols: Optical depth at reference wavelenght"
924	aeronlay_tauref=1.0E-1
925	call getin_p("aeronlay_tauref",aeronlay_tauref)
926	if (is_master) write(,)trim(rname)//&
927	": aeronlay_tauref = ",aeronlay_tauref
928
929	if (is_master) write(,)trim(rname)//&
930	": Generic n-layer aerosols: Reference wavelenght for optical depths (m)"
931	aeronlay_lamref=0.6E-6
932	call getin_p("aeronlay_lamref",aeronlay_lamref)
933	if (is_master) write(,)trim(rname)//&
934	": aeronlay_lamref = ",aeronlay_lamref
935
936	if (is_master) then
937	write(,)trim(rname)//&
938	": Generic n-layer aerosols: Vertical profile choice : "
939	write(,)trim(rname)//&
940	" (1) Tau btwn ptop and pbot follows atm. scale height"
941	write(,)trim(rname)//&
942	" or (2) Tau above pbot follows its own scale height"
943	endif
944	aeronlay_choice=1
945	call getin_p("aeronlay_choice",aeronlay_choice)
946	if (is_master) write(,)trim(rname)//&
947	": aeronlay_choice = ",aeronlay_choice
948
949	if (is_master) write(,)trim(rname)//&
950	": Generic n-layer aerosols: bottom pressures (Pa)"
951	aeronlay_pbot=2000.0
952	call getin_p("aeronlay_pbot",aeronlay_pbot)
953	if (is_master) write(,)trim(rname)//": aeronlay_pbot = ",aeronlay_pbot
954
955	if (is_master) write(,)trim(rname)//&
956	": Generic n-layer aerosols: (if choice=1) Top pressures (Pa) "
957	aeronlay_ptop=300000.0
958	call getin_p("aeronlay_ptop",aeronlay_ptop)
959	if (is_master) write(,)trim(rname)//": aeronlay_ptop = ",aeronlay_ptop
960
961	if (is_master) write(,)trim(rname)//&
962	": Generic n-layer aerosols: (if choice=2) Scale height / atm. scale height"
963	aeronlay_sclhght=0.2
964	call getin_p("aeronlay_sclhght",aeronlay_sclhght)
965	if (is_master) write(,)trim(rname)//&
966	": aeronlay_sclhght = ",aeronlay_sclhght
967
968	if (is_master) write(,)trim(rname)//&
969	": Generic n-layer aerosols: particles effective radii (m)"
970	aeronlay_size=1.e-6
971	call getin_p("aeronlay_size",aeronlay_size)
972	if (is_master) write(,)trim(rname)//": aeronlay_size = ",aeronlay_size
973
974	if (is_master) write(,)trim(rname)//&
975	": Generic n-layer aerosols: particles radii effective variance"
976	aeronlay_nueff=0.1
977	call getin_p("aeronlay_nueff",aeronlay_nueff)
978	if (is_master) write(,)trim(rname)//": aeronlay_nueff = ",aeronlay_nueff
979
980	if (is_master) write(,)trim(rname)//&
981	": Generic n-layer aerosols: VIS optical properties file"
982	optprop_aeronlay_vis = 'optprop_saturn_vis_n20.dat'
983	call getin_p("optprop_aeronlay_vis",optprop_aeronlay_vis)
984	if (is_master) write(,)trim(rname)//&
985	": optprop_aeronlay_vis = ",optprop_aeronlay_vis
986
987	if (is_master) write(,)trim(rname)//&
988	": Generic n-layer aerosols: IR optical properties file"
989	optprop_aeronlay_ir = 'optprop_saturn_ir_n20.dat'
990	call getin_p("optprop_aeronlay_ir",optprop_aeronlay_ir)
991	if (is_master) write(,)trim(rname)//&
992	": optprop_aeronlay_ir = ",optprop_aeronlay_ir
993
994
995	!=================================
996
997	if (is_master) write(,)trim(rname)//&
998	": Cloud pressure level (with kastprof only):"
999	cloudlvl=0. ! default value
1000	call getin_p("cloudlvl",cloudlvl)
1001	if (is_master) write(,)trim(rname)//": cloudlvl = ",cloudlvl
1002
1003	if (is_master) write(,)trim(rname)//&
1004	": Is the variable gas species radiatively active?"
1005	Tstrat=167.0
1006	varactive=.false.
1007	call getin_p("varactive",varactive)
1008	if (is_master) write(,)trim(rname)//": varactive = ",varactive
1009
1010	if (is_master) write(,)trim(rname)//&
1011	": Is the variable gas species distribution set?"
1012	varfixed=.false.
1013	call getin_p("varfixed",varfixed)
1014	if (is_master) write(,)trim(rname)//": varfixed = ",varfixed
1015
1016	if (is_master) write(,)trim(rname)//&
1017	": What is the saturation % of the variable species?"
1018	satval=0.8
1019	call getin_p("satval",satval)
1020	if (is_master) write(,)trim(rname)//": satval = ",satval
1021
1022
1023	! Test of incompatibility:
1024	! if varactive, then varfixed should be false
1025	if (varactive.and.varfixed) then
1026	call abort_physic(rname,'if varactive, varfixed must be OFF!',1)
1027	endif
1028
1029	if (is_master) write(,)trim(rname)//": Gravitationnal sedimentation ?"
1030	sedimentation=.false. ! default value
1031	call getin_p("sedimentation",sedimentation)
1032	if (is_master) write(,)trim(rname)//": sedimentation = ",sedimentation
1033
1034	if (is_master) write(,)trim(rname)//": Generic Condensation of tracers ?"
1035	generic_condensation=.false. !default value
1036	call getin_p("generic_condensation",generic_condensation)
1037	if (is_master) write(,)trim(rname)//": generic_condensation = ",generic_condensation
1038
1039	if (is_master) write(,)trim(rname)//": Generic rain of tracers ?"
1040	generic_rain=.false. !default value
1041	call getin_p("generic_rain",generic_rain)
1042	if (is_master) write(,)trim(rname)//": generic_rain = ",generic_rain
1043
1044	if (is_master) write(,)trim(rname)//": Include moist adjustement for GCS ?"
1045	moistadjustment_generic=.false. ! default value
1046	call getin_p("moistadjustment_generic",moistadjustment_generic)
1047	if (is_master) write(,)trim(rname)//": moistadjustment_generic = ", moistadjustment_generic
1048
1049	if (is_master) write(,)trim(rname)//": Moist convection inhibition for GCS ?"
1050	moist_convection_inhibition=.false. ! default value
1051	call getin_p("moist_convection_inhibition",moist_convection_inhibition)
1052	if (is_master) write(,)trim(rname)//": moist_convection_inhibition = ", moist_convection_inhibition
1053
1054	if (is_master) write(,)trim(rname)//": Virtual correction ?"
1055	virtual_correction=.false. !default value
1056	call getin_p("virtual_correction",virtual_correction)
1057	if (is_master) write(,)trim(rname)//": virtual_correction = ",virtual_correction
1058
1059	if (is_master) write(,)trim(rname)//": Compute water cycle ?"
1060	water=.false. ! default value
1061	call getin_p("water",water)
1062	if (is_master) write(,)trim(rname)//": water = ",water
1063
1064	! Test of incompatibility:
1065	! if water is true, there should be at least a tracer
1066	if (water.and.(.not.tracer)) then
1067	call abort_physic(rname,'if water is ON, tracer must be ON too!',1)
1068	endif
1069
1070	if (is_master) write(,)trim(rname)//": Include water condensation ?"
1071	watercond=.false. ! default value
1072	call getin_p("watercond",watercond)
1073	if (is_master) write(,)trim(rname)//": watercond = ",watercond
1074
1075	! Test of incompatibility:
1076	! if watercond is used, then water should be used too
1077	if (watercond.and.(.not.water)) then
1078	call abort_physic(rname,&
1079	'if watercond is used, water should be used too',1)
1080	endif
1081
1082	if (is_master) write(,)trim(rname)//": Include water moist adjustement ?"
1083	moistadjustment=.true. ! default value
1084	call getin_p("moistadjustment",moistadjustment)
1085	if (is_master) write(,)trim(rname)//": moistadjustment = ", moistadjustment
1086
1087	if (is_master) write(,)trim(rname)//": Include water precipitation ?"
1088	waterrain=.false. ! default value
1089	call getin_p("waterrain",waterrain)
1090	if (is_master) write(,)trim(rname)//": waterrain = ",waterrain
1091
1092	if (is_master) write(,)trim(rname)//": Include surface hydrology ?"
1093	hydrology=.false. ! default value
1094	call getin_p("hydrology",hydrology)
1095	if (is_master) write(,)trim(rname)//": hydrology = ",hydrology
1096
1097	if (is_master) write(,)trim(rname)//": Spectral Dependant albedo ?"
1098	albedo_spectral_mode=.false. ! default value
1099	call getin_p("albedo_spectral_mode",albedo_spectral_mode)
1100	if (is_master) write(,)trim(rname)//&
1101	": albedo_spectral_mode = ",albedo_spectral_mode
1102
1103	if (is_master) then
1104	write(,)trim(rname)//": Snow albedo ?"
1105	write(,)trim(rname)//": If albedo_spectral_mode=.true., then this "
1106	write(,)trim(rname)//": corresponds to the 0.5 microns snow albedo."
1107	endif
1108	albedosnow=0.5 ! default value
1109	call getin_p("albedosnow",albedosnow)
1110	if (is_master) write(,)trim(rname)//": albedosnow = ",albedosnow
1111
1112	if (is_master) write(,)trim(rname)//": Ocean albedo ?"
1113	alb_ocean=0.07 ! default value
1114	call getin_p("alb_ocean",alb_ocean)
1115	if (is_master) write(,)trim(rname)//": alb_ocean = ",alb_ocean
1116
1117	if (is_master) write(,)trim(rname)//": CO2 ice albedo ?"
1118	albedoco2ice=0.5 ! default value
1119	call getin_p("albedoco2ice",albedoco2ice)
1120	if (is_master) write(,)trim(rname)//": albedoco2ice = ",albedoco2ice
1121
1122	if (is_master) write(,)trim(rname)//": Maximum ice thickness ?"
1123	maxicethick=2.0 ! default value
1124	call getin_p("maxicethick",maxicethick)
1125	if (is_master) write(,)trim(rname)//": maxicethick = ",maxicethick
1126
1127	if (is_master) write(,)trim(rname)//": Freezing point of seawater ?"
1128	Tsaldiff=-1.8 ! default value
1129	call getin_p("Tsaldiff",Tsaldiff)
1130	if (is_master) write(,)trim(rname)//": Tsaldiff = ",Tsaldiff
1131
1132	if (is_master) write(,)trim(rname)//": Minimum eddy mix coeff in 1D ?"
1133	kmixmin=1.0e-2 ! default value
1134	call getin_p("kmixmin",kmixmin)
1135	if (is_master) write(,)trim(rname)//": kmixmin = ",kmixmin
1136
1137	varspec=.false. !default value
1138	call getin_p("varspec",varspec)
1139	if (varspec) then
1140	call getin_p("varspec_data",varspec_data)
1141	call getin_p("nvarlayer",nvarlayer)
1142	endif
1143
1144	if (is_master) write(,)'Predefined Cp from dynamics is ',cpp,'J kg^-1 K^-1'
1145	if (is_master) write(,)'Predefined Mg from dynamics is ',mugaz,'amu'
1146
1147	force_cpp=.false. ! default value
1148	call getin_p("force_cpp",force_cpp)
1149	if (force_cpp) then
1150	if (is_master) write(,)trim(rname)//": force_cpp = ",force_cpp
1151	if (is_master) write(,)trim(rname)//": force_cpp is deprecated.",&
1152	"Set cpp_mugaz_mode=1 in callfis to emulate force_cpp=.true."
1153	call abort_physic(rname,"Anyway, you need to set force_cpp=.false. to continue.",1)
1154	endif
1155
1156	if (is_master) write(,)trim(rname)//&
1157	": where do you want your cpp/mugaz value to come from?",&
1158	"=> 0: dynamics (3d), 1: forced in callfis (1d), 2: computed from gases.def (1d)?"
1159	cpp_mugaz_mode = 0 ! default value
1160	call getin_p("cpp_mugaz_mode",cpp_mugaz_mode)
1161	if (is_master) write(,)trim(rname)//": cpp_mugaz_mode = ",cpp_mugaz_mode
1162
1163	if ((cpp_mugaz_mode >= 1).and.(is_master).and.(ngrid>1)) then
1164	write(,)' !!!!! Be aware that having different values for cpp and mugaz in the dynamics and physics'
1165	write(,)' in 3D can result in very pathological behavior. You have been warned !!!!!'
1166	else if ((cpp_mugaz_mode == 0).and.(is_master).and.(ngrid==1)) then
1167	! for this specific 1D error we will remove run.def before aborting JL22
1168	call system("rm -rf run.def")
1169	call abort_physic(rname,"cpp_mugaz_mode must be >= 1 in 1d",1)
1170	endif
1171
1172	if (cpp_mugaz_mode == 1) then
1173	mugaz = -99999.
1174	if (is_master) write(,)trim(rname)//&
1175	": MEAN MOLECULAR MASS in g mol-1 ?"
1176	call getin_p("mugaz",mugaz)
1177	IF (mugaz.eq.-99999.) THEN
1178	call abort_physic(rname,"mugaz must be set if cpp_mugaz_mode = 1",1)
1179	ENDIF
1180	cpp = -99999.
1181	if (is_master) write(,)trim(rname)//&
1182	": SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
1183	call getin_p("cpp",cpp)
1184	IF (cpp.eq.-99999.) THEN
1185	PRINT *, "cpp must be set if cpp_mugaz_mode = 1"
1186	STOP
1187	ENDIF
1188	if (is_master) write(,)'New Cp from callfis is ',cpp,'J kg^-1 K^-1'
1189	if (is_master) write(,)'New Mg from callfis is ',mugaz,'amu'
1190
1191	endif ! of if (cpp_mugaz_mode == 1)
1192	call su_gases
1193	call calc_cpp_mugaz
1194
1195	if (is_master) then
1196	PRINT*,'--------------------------------------------'
1197	PRINT*
1198	PRINT*
1199	endif
1200	ELSE
1201	call abort_physic(rname,'Cannot read file callphys.def. Is it here ?',1)
1202	ENDIF ! of IF(iscallphys)
1203
1204	if (is_master) then
1205	PRINT*
1206	PRINT*,'inifis: daysec',daysec
1207	PRINT*
1208	PRINT*,'inifis: The radiative transfer is computed:'
1209	PRINT*,' each ',iradia,' physical time-step'
1210	PRINT,' or each ',iradiadtphys,' seconds'
1211	PRINT*
1212	endif
1213
1214	!-----------------------------------------------------------------------
1215	! Some more initialization:
1216	! ------------------------
1217
1218	! Initializations for comgeomfi_h
1219	#ifndef MESOSCALE
1220	totarea=SSUM(ngrid,parea,1)
1221	call planetwide_sumval(parea,totarea_planet)
1222
1223	!! those are defined in comdiurn_h.F90
1224	IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
1225	IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
1226	IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
1227	IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
1228
1229	DO ig=1,ngrid
1230	sinlat(ig)=sin(plat(ig))
1231	coslat(ig)=cos(plat(ig))
1232	sinlon(ig)=sin(plon(ig))
1233	coslon(ig)=cos(plon(ig))
1234	ENDDO
1235	#endif
1236	! initialize variables in radinc_h
1237	call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
1238
1239	! initialize variables and allocate arrays in radcommon_h
1240	call ini_radcommon_h(naerkind)
1241
1242	! allocate "comsoil_h" arrays
1243	call ini_comsoil_h(ngrid)
1244
1245	END SUBROUTINE inifis
1246
1247	END MODULE inifis_mod

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

Download in other formats:

Original Format