Context Navigation

inifis_mod.F90 @ 3567

Last change on this file since 3567 was 3437, checked in by emillour, 3 months ago

Generic PCM Integrate update from newton branch by LT:

remove call to old newtrelax routine and replace it with a new one, newtcool_MOCHA, tuned to the MOCHA intercomparison project

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

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

Download in other formats:

Original Format