Context Navigation

inifis_mod.F90 @ 3662

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

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

Download in other formats:

Original Format