Context Navigation

inifis_mod.F90 @ 3716

Last change on this file since 3716 was 3701, checked in by mmaurice, 3 months ago

Generic PCM

Change default reactfile from chemnetwork/reactfile to network.def.
You can still specify it callphys.def using the flag reactfile.

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

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

Download in other formats:

Original Format