Context Navigation

inifis_mod.F90 @ 3558

Last change on this file since 3558 was 3557, checked in by debatzbr, 4 weeks ago
Miscellaneous cleans + Set-up the physics for the implementation of the microphysical model.
File size: 59.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_n2
15	use datafile_mod, only: datadir,config_mufi,hazeprop_file,hazerad_file,hazemmr_file,hazedens_file
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 planetwide_mod, only: planetwide_sumval
24	use callkeys_mod
25	use surfdat_h
26	use wstats_mod, only: callstats
27	use ioipsl_getin_p_mod, only : getin_p
28	use mod_phys_lmdz_para, only : is_parallel, is_master, bcast
29
30	!=======================================================================
31	!
32	! purpose:
33	! -------
34	!
35	! Initialisation for the physical parametrisations of the LMD
36	! Generic Model.
37	!
38	! author: Frederic Hourdin 15 / 10 /93
39	! -------
40	! modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
41	! Ehouarn Millour (oct. 2008) tracers are now identified
42	! by their names and may not be contiguously
43	! stored in the q(:,:,:,:) array
44	! E.M. (june 2009) use getin routine to load parameters
45	!
46	!
47	! arguments:
48	! ----------
49	!
50	! input:
51	! ------
52	!
53	! ngrid Size of the horizontal grid.
54	! All internal loops are performed on that grid.
55	! nlayer Number of vertical layers.
56	! pdayref Day of reference for the simulation
57	! pday Number of days counted from the North. Spring
58	! equinoxe.
59	!
60	!=======================================================================
61	!
62	!-----------------------------------------------------------------------
63	! declarations:
64	! -------------
65	use datafile_mod, only: datadir
66	use ioipsl_getin_p_mod, only: getin_p
67	IMPLICIT NONE
68
69
70
71	REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep
72	INTEGER,INTENT(IN) :: nday
73	INTEGER,INTENT(IN) :: ngrid,nlayer,nq
74	REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid)
75	integer,intent(in) :: day_ini
76
77	INTEGER ig
78	CHARACTER(len=20) :: rname="inifis" ! routine name, for messages
79
80	EXTERNAL iniorbit,orbite
81	EXTERNAL SSUM
82	REAL SSUM
83
84	! Initialize flags lunout, prt_level, debug (in print_control_mod)
85	CALL init_print_control
86
87	! initialize constants in comcstfi_mod
88	rad=prad
89	cpp=pcpp
90	g=pg
91	r=pr
92	rcp=r/cpp
93	mugaz=8.314*1000./pr ! dummy init
94	pi=2.*asin(1.)
95	avocado = 6.02214179e23 ! added by RW
96
97	! Initialize some "temporal and calendar" related variables
98	#ifndef MESOSCALE
99	CALL init_time(day_ini,pdaysec,nday,ptimestep)
100	#endif
101
102	! read in some parameters from "run.def" for physics,
103	! or shared between dynamics and physics.
104	call getin_p("ecritphy",ecritphy) ! frequency of outputs in physics,
105	! in dynamical steps
106	call getin_p("day_step",day_step) ! number of dynamical steps per day
107	call getin_p("iphysiq",iphysiq) ! call physics every iphysiq dyn step
108
109	! do we read a startphy.nc file? (default: .true.)
110	call getin_p("startphy_file",startphy_file)
111
112	! --------------------------------------------------------------
113	! Reading the "callphys.def" file controlling some key options
114	! --------------------------------------------------------------
115
116	IF (is_master) THEN
117	! check if 'callphys.def' file is around
118	inquire(file="callphys.def",exist=iscallphys)
119	write(,) trim(rname)//": iscallphys=",iscallphys
120	ENDIF
121	call bcast(iscallphys)
122
123	IF(iscallphys) THEN
124
125	if (is_master) then
126	write(,)
127	write(,)'--------------------------------------------'
128	write(,)trim(rname)//': Parameters for the physics (callphys.def)'
129	write(,)'--------------------------------------------'
130	endif
131
132	if (is_master) write(,) trim(rname)//&
133	": Directory where external input files are:"
134	! default 'datadir' is set in "datadir_mod"
135	call getin_p("datadir",datadir) ! default path
136	if (is_master) write(,) trim(rname)//": datadir = ",trim(datadir)
137
138	if (is_master) write(,) "Haze optical properties datafile"
139	hazeprop_file="optprop_tholins_fractal050nm.dat" ! default file
140	call getin_p("hazeprop_file",hazeprop_file)
141	if (is_master) write(,) trim(rname)//" hazeprop_file = ",trim(hazeprop_file)
142
143	if (is_master) write(,) "Haze radii datafile"
144	hazerad_file="hazerad.txt" ! default file
145	call getin_p("hazerad_file",hazerad_file)
146	if (is_master) write(,) trim(rname)//" hazerad_file = ",trim(hazerad_file)
147
148	if (is_master) write(,) "Haze mmr datafile"
149	hazemmr_file="None" ! default file
150	call getin_p("hazemmr_file",hazemmr_file)
151	if (is_master) write(,) trim(rname)//" hazemmr_file = ",trim(hazemmr_file)
152	if (is_master) write(,) "Haze density datafile"
153	hazedens_file="None" ! default file
154	call getin_p("hazedens_file",hazedens_file)
155	if (is_master) write(,) trim(rname)//" hazedens_file = ",trim(hazedens_file)
156
157	if (is_master) write(,) trim(rname)//&
158	": Run with or without tracer transport ?"
159	tracer=.false. ! default value
160	call getin_p("tracer",tracer)
161	if (is_master) write(,) trim(rname)//": tracer = ",tracer
162
163	if (is_master) write(,) trim(rname)//&
164	": Run with or without atm mass update "//&
165	" due to tracer evaporation/condensation?"
166	mass_redistrib=.false. ! default value
167	call getin_p("mass_redistrib",mass_redistrib)
168	if (is_master) write(,) trim(rname)//&
169	": mass_redistrib = ",mass_redistrib
170
171	if (is_master) then
172	write(,) trim(rname)//": Diurnal cycle ?"
173	write(,) trim(rname)//&
174	": (if diurnal=false, diurnal averaged solar heating)"
175	endif
176	diurnal=.true. ! default value
177	call getin_p("diurnal",diurnal)
178	if (is_master) write(,) trim(rname)//": diurnal = ",diurnal
179
180	if (is_master) then
181	write(,) trim(rname)//": Seasonal cycle ?"
182	write(,) trim(rname)//&
183	": (if season=false, Ls stays constant, to value ", &
184	"set in 'start'"
185	endif
186	season=.true. ! default value
187	call getin_p("season",season)
188	if (is_master) write(,) trim(rname)//": season = ",season
189
190	if (is_master) then
191	write(,) trim(rname)//": Fast mode (nogcm) ?"
192	endif
193	fast=.false. ! default value
194	call getin_p("fast",fast)
195	if (is_master) write(,) trim(rname)//": fast = ",fast
196
197	if (is_master) write(,) trim(rname)//"Run with Triton orbit ?"
198	triton=.false. ! default value
199	call getin_p("triton",triton)
200	if (is_master) write(,) trim(rname)//" triton = ",triton
201	convergeps=.false. ! default value
202	call getin_p("convergeps",convergeps)
203	if (is_master) write(,) trim(rname)//" convergeps = ",convergeps
204
205	conservn2=.false. ! default value
206	call getin_p("conservn2",conservn2)
207	if (is_master) write(,) trim(rname)//" conservn2 = ",conservn2
208	conservch4=.false. ! default value
209	call getin_p("conservch4",conservch4)
210	if (is_master) write(,) trim(rname)//" conservch4 = ",conservch4
211
212	if (is_master) write(,) trim(rname)//"KBO runs (eris, makemake) ?"
213	kbo=.false. ! default value
214	call getin_p("kbo",kbo)
215	if (is_master) write(,) trim(rname)//" kbo = ",kbo
216
217	if (is_master) write(,) trim(rname)//"Specific paleo run ?"
218	paleo=.false. ! default value
219	call getin_p("paleo",paleo)
220	if (is_master) write(,) trim(rname)//" paleo = ",paleo
221
222	if (is_master) write(,) trim(rname)//"paleoclimate step (Earth years) "
223	paleoyears=10000. ! default value
224	call getin_p("paleoyears",paleoyears)
225	if (is_master) write(,) trim(rname)//" paleoyears = ",paleoyears
226
227	if (is_master) write(,) trim(rname)//&
228	": No seasonal cycle: initial day to lock the run during restart"
229	noseason_day=0.0 ! default value
230	call getin_p("noseason_day",noseason_day)
231	if (is_master) write(,) trim(rname)//": noseason_day=", noseason_day
232
233	if (is_master) write(,) trim(rname)//&
234	": Write some extra output to the screen ?"
235	lwrite=.false. ! default value
236	call getin_p("lwrite",lwrite)
237	if (is_master) write(,) trim(rname)//": lwrite = ",lwrite
238
239	if (is_master) write(,) trim(rname)//&
240	": Save statistics in file stats.nc ?"
241	callstats=.true. ! default value
242	call getin_p("callstats",callstats)
243	if (is_master) write(,) trim(rname)//": callstats = ",callstats
244
245	if (is_master) write(,) trim(rname)//&
246	": Test energy conservation of model physics ?"
247	enertest=.false. ! default value
248	call getin_p("enertest",enertest)
249	if (is_master) write(,) trim(rname)//": enertest = ",enertest
250
251	if (is_master) write(,) trim(rname)//": call radiative transfer ?"
252	callrad=.true. ! default value
253	call getin_p("callrad",callrad)
254	if (is_master) write(,) trim(rname)//": callrad = ",callrad
255
256	if (is_master) write(,) trim(rname)//&
257	": call correlated-k radiative transfer ?"
258	corrk=.true. ! default value
259	call getin_p("corrk",corrk)
260	if (is_master) write(,) trim(rname)//": corrk = ",corrk
261
262	if (is_master) write(,) trim(rname)//&
263	": prohibit calculations outside corrk T grid?"
264	strictboundcorrk=.true. ! default value
265	call getin_p("strictboundcorrk",strictboundcorrk)
266	if (is_master) write(,) trim(rname)//&
267	": strictboundcorrk = ",strictboundcorrk
268
269	if (is_master) write(,) trim(rname)//&
270	": prohibit calculations outside CIA T grid?"
271	strictboundcia=.true. ! default value
272	call getin_p("strictboundcia",strictboundcia)
273	if (is_master) write(,) trim(rname)//&
274	": strictboundcia = ",strictboundcia
275
276	if (is_master) write(,) trim(rname)//&
277	": Minimum atmospheric temperature for Planck function integration ?"
278	tplanckmin=30.0 ! default value
279	call getin_p("tplanckmin",tplanckmin)
280	if (is_master) write(,) trim(rname)//": tplanckmin = ",tplanckmin
281
282	if (is_master) write(,) trim(rname)//&
283	": Maximum atmospheric temperature for Planck function integration ?"
284	tplanckmax=1500.0 ! default value
285	call getin_p("tplanckmax",tplanckmax)
286	if (is_master) write(,) trim(rname)//": tplanckmax = ",tplanckmax
287
288	if (is_master) write(,) trim(rname)//&
289	": Temperature step for Planck function integration ?"
290	dtplanck=0.1 ! default value
291	call getin_p("dtplanck",dtplanck)
292	if (is_master) write(,) trim(rname)//": dtplanck = ",dtplanck
293
294	if (is_master) write(,) trim(rname)//&
295	": call gaseous absorption in the visible bands?"//&
296	" (matters only if callrad=T)"
297	callgasvis=.false. ! default value
298	call getin_p("callgasvis",callgasvis)
299	if (is_master) write(,) trim(rname)//": callgasvis = ",callgasvis
300
301	if (is_master) write(,) trim(rname)//&
302	": call continuum opacities in radiative transfer ?"//&
303	" (matters only if callrad=T)"
304	continuum=.true. ! default value
305	call getin_p("continuum",continuum)
306	if (is_master) write(,) trim(rname)//": continuum = ",continuum
307
308	if (is_master) write(,) trim(rname)//&
309	": call turbulent vertical diffusion ?"
310	calldifv=.false. ! default value
311	call getin_p("calldifv",calldifv)
312	if (is_master) write(,) trim(rname)//": calldifv = ",calldifv
313
314	if (is_master) write(,) trim(rname)//": use turbdiff instead of vdifc ?"
315	UseTurbDiff=.false. ! default value
316	call getin_p("UseTurbDiff",UseTurbDiff)
317	if (is_master) write(,) trim(rname)//": UseTurbDiff = ",UseTurbDiff
318
319	if (is_master) write(,) trim(rname)//": call convective adjustment ?"
320	calladj=.false. ! default value
321	call getin_p("calladj",calladj)
322	if (is_master) write(,) trim(rname)//": calladj = ",calladj
323
324	if (is_master) write(,) trim(rname)//&
325	": Radiative timescale for Newtonian cooling (in Earth days)?"
326	tau_relax=30. ! default value
327	call getin_p("tau_relax",tau_relax)
328	if (is_master) write(,) trim(rname)//": tau_relax = ",tau_relax
329	tau_relax=tau_relax243600 ! convert Earth days --> seconds
330
331	if (is_master) write(,)trim(rname)//&
332	": call thermal conduction in the soil ?"
333	callsoil=.true. ! default value
334	call getin_p("callsoil",callsoil)
335	if (is_master) write(,) trim(rname)//": callsoil = ",callsoil
336
337	if (is_master) write(,)trim(rname)//&
338	": Rad transfer is computed every iradia", &
339	" physical timestep"
340	iradia=1 ! default value
341	call getin_p("iradia",iradia)
342	if (is_master) write(,)trim(rname)//": iradia = ",iradia
343
344	if (is_master) write(,)trim(rname)//": Rayleigh scattering ?"
345	rayleigh=.false.
346	call getin_p("rayleigh",rayleigh)
347	if (is_master) write(,)trim(rname)//": rayleigh = ",rayleigh
348
349	if (is_master) write(,) trim(rname)//&
350	": Use blackbody for stellar spectrum ?"
351	stelbbody=.false. ! default value
352	call getin_p("stelbbody",stelbbody)
353	if (is_master) write(,) trim(rname)//": stelbbody = ",stelbbody
354
355	if (is_master) write(,) trim(rname)//": Stellar blackbody temperature ?"
356	stelTbb=5800.0 ! default value
357	call getin_p("stelTbb",stelTbb)
358	if (is_master) write(,) trim(rname)//": stelTbb = ",stelTbb
359
360	if (is_master) write(,) trim(rname)//": Output mean OLR in 1D?"
361	meanOLR=.false.
362	call getin_p("meanOLR",meanOLR)
363	if (is_master) write(,) trim(rname)//": meanOLR = ",meanOLR
364
365	if (is_master) write(,)trim(rname)//": Output spectral OLR in 3D?"
366	specOLR=.false.
367	call getin_p("specOLR",specOLR)
368	if (is_master) write(,)trim(rname)//": specOLR = ",specOLR
369
370	if (is_master) write(,)trim(rname)//&
371	": Output diagnostic optical thickness attenuation exp(-tau)?"
372	diagdtau=.false.
373	call getin_p("diagdtau",diagdtau)
374	if (is_master) write(,)trim(rname)//": diagdtau = ",diagdtau
375
376	if (is_master) write(,)trim(rname)//": Operate in kastprof mode?"
377	kastprof=.false.
378	call getin_p("kastprof",kastprof)
379	if (is_master) write(,)trim(rname)//": kastprof = ",kastprof
380
381	if (is_master) write(,)trim(rname)//&
382	": Uniform absorption in radiative transfer?"
383	graybody=.false.
384	call getin_p("graybody",graybody)
385	if (is_master) write(,)trim(rname)//": graybody = ",graybody
386
387	! Soil model
388	if (is_master) write(,)trim(rname)//&
389	": Number of sub-surface layers for soil scheme?"
390	! default value of nsoilmx set in comsoil_h
391	call getin_p("nsoilmx",nsoilmx)
392	if (is_master) write(,)trim(rname)//": nsoilmx=",nsoilmx
393
394	if (is_master) write(,)trim(rname)//&
395	": Thickness of topmost soil layer (m)?"
396	! default value of lay1_soil set in comsoil_h
397	call getin_p("lay1_soil",lay1_soil)
398	if (is_master) write(,)trim(rname)//": lay1_soil = ",lay1_soil
399
400	if (is_master) write(,)trim(rname)//&
401	": Coefficient for soil layer thickness distribution?"
402	! default value of alpha_soil set in comsoil_h
403	call getin_p("alpha_soil",alpha_soil)
404	if (is_master) write(,)trim(rname)//": alpha_soil = ",alpha_soil
405
406	! Chemistry in the thermosphere
407	if (is_master) write(,) trim(rname)//": Use deposition ?"
408	depos=.false. ! default value
409	call getin_p("depos",depos)
410	if (is_master) write(,) trim(rname)//": depos = ",depos
411
412	if (is_master) write(,)trim(rname)//": Production of haze ?"
413	haze=.false. ! default value
414	call getin_p("haze",haze)
415	if (is_master) write(,)trim(rname)//": haze = ",haze
416
417
418
419	if (is_master) write(,)trim(rname)// "call thermal conduction ?"
420	callconduct=.false. ! default value
421	call getin_p("callconduct",callconduct)
422	if (is_master) write(,)trim(rname)// " callconduct = ",callconduct
423
424	if (is_master) write(,)trim(rname)// "call phitop ?"
425	phitop=0. ! default value
426	call getin_p("phitop",phitop)
427	if (is_master) write(,)trim(rname)// " phitop = ",phitop
428
429	if (is_master) write(,)trim(rname)// "call molecular viscosity ?"
430	callmolvis=.false. ! default value
431	call getin_p("callmolvis",callmolvis)
432	if (is_master) write(,)trim(rname)// " callmolvis = ",callmolvis
433
434	if (is_master) write(,)trim(rname)// "call molecular diffusion ?"
435	callmoldiff=.false. ! default value
436	call getin_p("callmoldiff",callmoldiff)
437	if (is_master) write(,)trim(rname)// " callmoldiff = ",callmoldiff
438
439	! Global1D mean and solar zenith angle
440
441	if (ngrid.eq.1) then
442	PRINT*, 'Simulate global averaged conditions ?'
443	global1d = .false. ! default value
444	call getin_p("global1d",global1d)
445	write(,) "global1d = ",global1d
446
447	! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
448	if (global1d.and.diurnal) then
449	call abort_physic(rname,'if global1d is true, diurnal must be set to false',1)
450	endif
451
452	if (global1d) then
453	PRINT *,'Solar Zenith angle (deg.) ?'
454	PRINT *,'(assumed for averaged solar flux S/4)'
455	szangle=60.0 ! default value
456	call getin_p("szangle",szangle)
457	write(,) "szangle = ",szangle
458	endif
459	endif ! of if (ngrid.eq.1)
460
461	! Test of incompatibility:
462	! if kastprof used, we must be in 1D
463	if (kastprof.and.(ngrid.gt.1)) then
464	call abort_physic(rname,'kastprof can only be used in 1D!',1)
465	endif
466
467	if (is_master) write(,)trim(rname)//&
468	": Stratospheric temperature for kastprof mode?"
469	Tstrat=167.0
470	call getin_p("Tstrat",Tstrat)
471	if (is_master) write(,)trim(rname)//": Tstrat = ",Tstrat
472
473	if (is_master) write(,)trim(rname)//": Remove lower boundary?"
474	nosurf=.false.
475	call getin_p("nosurf",nosurf)
476	if (is_master) write(,)trim(rname)//": nosurf = ",nosurf
477
478	! Tests of incompatibility:
479	if (nosurf.and.callsoil) then
480	if (is_master) then
481	write(,)trim(rname)//'nosurf not compatible with soil scheme!'
482	write(,)trim(rname)//'... got to make a choice!'
483	endif
484	call abort_physic(rname,"nosurf not compatible with soil scheme!",1)
485	endif
486
487	if (is_master) write(,)trim(rname)//": Add an internal heat flux?", &
488	"... matters only if callsoil=F"
489	intheat=0.
490	call getin_p("intheat",intheat)
491	if (is_master) write(,)trim(rname)//": intheat = ",intheat
492
493
494	if (is_master) write(,)trim(rname)//": Test physics timescale in 1D?"
495	testradtimes=.false.
496	call getin_p("testradtimes",testradtimes)
497	if (is_master) write(,)trim(rname)//": testradtimes = ",testradtimes
498
499	! Test of incompatibility:
500	! if testradtimes used, we must be in 1D
501	if (testradtimes.and.(ngrid.gt.1)) then
502	call abort_physic(rname,'testradtimes can only be used in 1D!',1)
503	endif
504
505	if (is_master) write(,)trim(rname)//": Default planetary temperature?"
506	tplanet=215.0
507	call getin_p("tplanet",tplanet)
508	if (is_master) write(,)trim(rname)//": tplanet = ",tplanet
509
510	if (is_master) write(,)trim(rname)//": Which star?"
511	startype=1 ! default value = Sol
512	call getin_p("startype",startype)
513	if (is_master) write(,)trim(rname)//": startype = ",startype
514
515	if (is_master) write(,)trim(rname)//": Value of stellar flux at 1 AU?"
516	Fat1AU=1356.0 ! default value = Sol today
517	call getin_p("Fat1AU",Fat1AU)
518	if (is_master) write(,)trim(rname)//": Fat1AU = ",Fat1AU
519
520
521	! TRACERS:
522
523	if (is_master) write(,)trim(rname)//": fixed radii for Cloud particles?"
524	radfixed=.false. ! default value
525	call getin_p("radfixed",radfixed)
526	if (is_master) write(,)trim(rname)//": radfixed = ",radfixed
527
528	if(kastprof)then
529	radfixed=.true.
530	endif
531
532	if (is_master) write(,)trim(rname)//&
533	"Number of radiatively active aerosols:"
534	naerkind=0 ! default value
535	call getin_p("naerkind",naerkind)
536	if (is_master) write(,)trim(rname)//": naerkind = ",naerkind
537
538
539	!***************************************************************
540	!! TRACERS options
541
542	if (is_master)write(,)trim(rname)//&
543	"call N2 condensation ?"
544	n2cond=.true. ! default value
545	call getin_p("n2cond",n2cond)
546	if (is_master)write(,)trim(rname)//&
547	" n2cond = ",n2cond
548
549	if (is_master)write(,)trim(rname)//&
550	"call N2 cloud condensation ?"
551	condensn2=.false. ! default value
552	call getin_p("condensn2",condensn2)
553	if (is_master)write(,)trim(rname)//&
554	"condensn2 = ",condensn2
555
556	if (is_master)write(,)trim(rname)//&
557	"call no N2 frost formation?"
558	no_n2frost=.false. ! default value
559	call getin_p("no_n2frost",no_n2frost)
560	if (is_master)write(,)trim(rname)//&
561	"no_n2frost = ",no_n2frost
562
563	if (is_master)write(,)trim(rname)//&
564	"N2 condensation subtimestep?"
565	nbsub=20 ! default value
566	call getin_p("nbsub",nbsub)
567	if (is_master)write(,)trim(rname)//&
568	" nbsub = ",nbsub
569
570	if (is_master)write(,)trim(rname)//&
571	"Gravitationnal sedimentation ?"
572	sedimentation=.true. ! default value
573	call getin_p("sedimentation",sedimentation)
574	if (is_master)write(,)trim(rname)//&
575	" sedimentation = ",sedimentation
576
577	if (is_master)write(,)trim(rname)//&
578	"Compute glacial flow ?"
579	glaflow=.false. ! default value
580	call getin_p("glaflow",glaflow)
581	if (is_master)write(,)trim(rname)//&
582	" glaflow = ",glaflow
583
584	if (is_master)write(,)trim(rname)//&
585	"Compute methane cycle ?"
586	methane=.false. ! default value
587	call getin_p("methane",methane)
588	if (is_master)write(,)trim(rname)//&
589	" methane = ",methane
590	condmetsurf=.true. ! default value
591	call getin_p("condmetsurf",condmetsurf)
592	if (is_master)write(,)trim(rname)//&
593	" condmetsurf = ",condmetsurf
594
595	if (is_master)write(,)trim(rname)//&
596	"Compute methane clouds ?"
597	metcloud=.false. ! default value
598	call getin_p("metcloud",metcloud)
599	if (is_master)write(,)trim(rname)//&
600	" metcloud = ",metcloud
601
602	if (is_master)write(,)trim(rname)//&
603	"Compute CO cycle ?"
604	carbox=.false. ! default value
605	call getin_p("carbox",carbox)
606	if (is_master)write(,)trim(rname)//&
607	" carbox = ",carbox
608	condcosurf=.true. ! default value
609	call getin_p("condcosurf",condcosurf)
610	if (is_master)write(,)trim(rname)//&
611	" condcosurf = ",condcosurf
612
613	if (is_master)write(,)trim(rname)//&
614	"Compute CO clouds ?"
615	monoxcloud=.false. ! default value
616	call getin_p("monoxcloud",monoxcloud)
617	if (is_master)write(,)trim(rname)//&
618	" monoxcloud = ",monoxcloud
619
620	if (is_master)write(,)trim(rname)//&
621	"atmospheric redistribution (s):"
622	tau_n2=1. ! default value
623	call getin_p("tau_n2",tau_n2)
624	if (is_master)write(,)trim(rname)//&
625	" tau_n2 = ",tau_n2
626	tau_ch4=1.E7 ! default value
627	call getin_p("tau_ch4",tau_ch4)
628	if (is_master)write(,)trim(rname)//&
629	" tau_ch4 = ",tau_ch4
630	tau_co=1. ! default value
631	call getin_p("tau_co",tau_co)
632	if (is_master)write(,)trim(rname)//&
633	" tau_co = ",tau_co
634	tau_prechaze=1. ! default value
635	call getin_p("tau_prechaze",tau_prechaze)
636	if (is_master)write(,)trim(rname)//&
637	" tau_prechaze = ",tau_prechaze
638
639	if (is_master)write(,)trim(rname)//&
640	"Day fraction for limited cold trap in SP?"
641	dayfrac=0. ! default value
642	call getin_p("dayfrac",dayfrac)
643	if (is_master)write(,)trim(rname)//&
644	" dayfrac = ",dayfrac
645	thresh_non2=0. ! default value
646	call getin_p("thresh_non2",thresh_non2)
647	if (is_master)write(,)trim(rname)//&
648	" thresh_non2 = ",thresh_non2
649
650	! use Pluto.old routines
651	if (is_master) write(,) trim(rname)//": use vdifc from old Pluto ?"
652	oldplutovdifc=.false. ! default value
653	call getin_p("oldplutovdifc",oldplutovdifc)
654	if (is_master) write(,) trim(rname)//": oldplutovdifc = ",oldplutovdifc
655
656	if (is_master) write(,) trim(rname)//&
657	": call pluto.old correlated-k radiative transfer ?"
658	oldplutocorrk=.false. ! default value
659	call getin_p("oldplutocorrk",oldplutocorrk)
660	if (is_master) write(,) trim(rname)//": oldplutocorrk = ",oldplutocorrk
661
662	if (is_master) write(,) trim(rname)//&
663	": call pluto.old sedimentation ?"
664	oldplutosedim=.false. ! default value
665	call getin_p("oldplutosedim",oldplutosedim)
666	if (is_master) write(,) trim(rname)//": oldplutosedim = ",oldplutosedim
667
668	!***************************************************************
669	!! Haze options
670
671	! Microphysical moment model
672	! ~~~~~~~~~~~~~~~~~~~~~~~~~~
673	if (is_master) write(,) "Run with or without microphysics?"
674	callmufi=.false. ! default value
675	call getin_p("callmufi",callmufi)
676	if (is_master) write(,)" callmufi = ",callmufi
677
678	! sanity check
679	if (callmufi.and.(.not.tracer)) then
680	print*,"You are running microphysics without tracer"
681	print*,"Please start again with tracer =.true."
682	stop
683	endif
684
685	if (is_master) write(,) "Path to microphysical config file?"
686	config_mufi='datagcm/microphysics/config.cfg' ! default value
687	call getin_p("config_mufi",config_mufi)
688	if (is_master) write(,)" config_mufi = ",config_mufi
689
690	if (is_master) write(,) "Use haze production from CH4 photolysis or production rate?"
691	call_haze_prod_pCH4=.false. ! default value
692	call getin_p("call_haze_prod_pCH4",call_haze_prod_pCH4)
693	if (is_master) write(,)" call_haze_prod_pCH4 = ",call_haze_prod_pCH4
694
695	if (is_master) write(,) "Pressure level of aerosols production (Pa)?"
696	haze_p_prod=1.0e-2 ! default value
697	call getin_p("haze_p_prod",haze_p_prod)
698	if (is_master) write(,)" haze_p_prod = ",haze_p_prod
699
700	if (is_master) write(,) "Aerosol production rate (kg.m-2.s-1)?"
701	haze_tx_prod=9.8e-14 ! default value
702	call getin_p("haze_tx_prod",haze_tx_prod)
703	if (is_master) write(,)" haze_tx_prod = ",haze_tx_prod
704
705	if (is_master) write(,) "Equivalent radius production (m)?"
706	haze_rc_prod=1.0e-9 ! default value
707	call getin_p("haze_rc_prod",haze_rc_prod)
708	if (is_master) write(,)" haze_rc_prod = ",haze_rc_prod
709
710	if (is_master) write(,) "Monomer radius (m)?"
711	haze_rm=1.0e-8 ! default value
712	call getin_p("haze_rm",haze_rm)
713	if (is_master) write(,)" haze_rm = ",haze_rm
714
715	if (is_master) write(,) "Aerosol's fractal dimension?"
716	haze_df=2.0 ! default value
717	call getin_p("haze_df",haze_df)
718	if (is_master) write(,)" haze_df = ",haze_df
719
720	if (is_master) write(,) "Aerosol density (kg.m-3)?"
721	haze_rho=800.0 ! default value
722	call getin_p("haze_rho",haze_rho)
723	if (is_master) write(,)" haze_rho = ",haze_rho
724
725	if (is_master) write(,) "Radius of air molecule (m)?"
726	air_rad=1.75e-10 ! default value
727	call getin_p("air_rad",air_rad)
728	if (is_master) write(,)" air_rad = ",air_rad
729
730	! Pluto haze model
731	! ~~~~~~~~~~~~~~~~
732	if (is_master)write(,)trim(rname)//&
733	"Production of haze ?"
734	haze=.false. ! default value
735	call getin_p("haze",haze)
736	if (is_master)write(,)trim(rname)//&
737	" haze = ",haze
738	hazeconservch4=.false. ! conservch4, by default value ch4 is photolyzed
739	call getin_p("hazeconservch4",hazeconservch4)
740	if (is_master)write(,)trim(rname)//&
741	"hazconservch4 = ",hazeconservch4
742	if (is_master)write(,)trim(rname)//&
743	"Production of haze (fast version) ?"
744	fasthaze=.false. ! default value
745	call getin_p("fasthaze",fasthaze)
746	if (is_master)write(,)trim(rname)//&
747	"fasthaze = ",fasthaze
748
749	if (is_master)write(,)trim(rname)//&
750	"Add source of haze ?"
751	source_haze=.false. ! default value
752	call getin_p("source_haze",source_haze)
753	if (is_master)write(,)trim(rname)//&
754	" source_haze = ",source_haze
755	mode_hs=0 ! mode haze source default value
756	call getin_p("mode_hs",mode_hs)
757	if (is_master)write(,)trim(rname)//&
758	" mode_hs",mode_hs
759	kfix=1 ! default value
760	call getin_p("kfix",kfix)
761	if (is_master)write(,)trim(rname)//&
762	"levels kfix",kfix
763	fracsource=0.1 ! default value
764	call getin_p("fracsource",fracsource)
765	if (is_master)write(,)trim(rname)//&
766	" fracsource",fracsource
767	latsource=30. ! default value
768	call getin_p("latsource",latsource)
769	if (is_master)write(,)trim(rname)//&
770	" latsource",latsource
771	lonsource=180. ! default value
772	call getin_p("lonsource",lonsource)
773	if (is_master)write(,)trim(rname)//&
774	" lonsource",lonsource
775
776	if (is_master)write(,)trim(rname)//&
777	"Radiatively active haze ?"
778	aerohaze=.false. ! default value
779	call getin_p("aerohaze",aerohaze)
780	if (is_master)write(,)trim(rname)//&
781	"aerohaze = ",aerohaze
782
783	if (is_master)write(,)trim(rname)//&
784	"Haze monomer radius ?"
785	rad_haze=10.e-9 ! default value
786	call getin_p("rad_haze",rad_haze)
787	if (is_master)write(,)trim(rname)//&
788	"rad_haze = ",rad_haze
789
790	if (is_master)write(,)trim(rname)//&
791	"fractal particle ?"
792	fractal=.false. ! default value
793	call getin_p("fractal",fractal)
794	if (is_master)write(,)trim(rname)//&
795	"fractal = ",fractal
796	nb_monomer=10 ! default value
797	call getin_p("nb_monomer",nb_monomer)
798	if (is_master)write(,)trim(rname)//&
799	"nb_monomer = ",nb_monomer
800
801	if (is_master)write(,)trim(rname)//&
802	"fixed gaseous CH4 mixing ratio for rad transfer ?"
803	ch4fix=.false. ! default value
804	call getin_p("ch4fix",ch4fix)
805	if (is_master)write(,)trim(rname)//&
806	" ch4fix = ",ch4fix
807	if (is_master)write(,)trim(rname)//&
808	"fixed gaseous CH4 mixing ratio for rad transfer ?"
809	vmrch4fix=0.5 ! default value
810	call getin_p("vmrch4fix",vmrch4fix)
811	if (is_master)write(,)trim(rname)//&
812	" vmrch4fix = ",vmrch4fix
813	vmrch4_proffix=.false. ! default value
814	call getin_p("vmrch4_proffix",vmrch4_proffix)
815	if (is_master)write(,)trim(rname)//&
816	" vmrch4_proffix = ",vmrch4_proffix
817
818	if (is_master)write(,)trim(rname)//&
819	"call specific cooling for radiative transfer ?"
820	cooling=.false. ! default value
821	call getin_p("cooling",cooling)
822	if (is_master)write(,)trim(rname)//&
823	" cooling = ",cooling
824
825	if (is_master)write(,)trim(rname)//&
826	"NLTE correction for methane heating rates?"
827	nlte=.false. ! default value
828	call getin_p("nlte",nlte)
829	if (is_master)write(,)trim(rname)//&
830	" nlte = ",nlte
831	strobel=.false. ! default value
832	call getin_p("strobel",strobel)
833	if (is_master)write(,)trim(rname)//&
834	" strobel = ",strobel
835
836	if (is_master)write(,)trim(rname)//&
837	"fixed radius profile from txt file ?"
838	haze_radproffix=.false. ! default value
839	call getin_p("haze_radproffix",haze_radproffix)
840	if (is_master)write(,)trim(rname)//&
841	"haze_radproffix = ",haze_radproffix
842	if (is_master)write(,)trim(rname)//&
843	"fixed MMR profile from txt file ?"
844	haze_proffix=.false. ! default value
845	call getin_p("haze_proffix",haze_proffix)
846	if (is_master)write(,)trim(rname)//&
847	"haze_proffix = ",haze_proffix
848
849	if (is_master)write(,)trim(rname)//&
850	"Number mix ratio of haze particles for co clouds:"
851	Nmix_co=100000. ! default value
852	call getin_p("Nmix_co",Nmix_co)
853	if (is_master)write(,)trim(rname)//&
854	" Nmix_co = ",Nmix_co
855
856	if (is_master)write(,)trim(rname)//&
857	"Number mix ratio of haze particles for ch4 clouds:"
858	Nmix_ch4=100000. ! default value
859	call getin_p("Nmix_ch4",Nmix_ch4)
860	if (is_master)write(,)trim(rname)//&
861	" Nmix_ch4 = ",Nmix_ch4
862
863	!***************************************************************
864	!! Surface properties
865
866	!********* N2 *******************************
867
868	if (is_master)write(,)trim(rname)//&
869	"Mode for changing N2 albedo"
870	mode_n2=0 ! default value
871	call getin_p("mode_n2",mode_n2)
872	if (is_master)write(,)trim(rname)//&
873	" mode_n2 = ",mode_n2
874	thres_n2ice=1. ! default value
875	call getin_p("thres_n2ice",thres_n2ice)
876	if (is_master)write(,)trim(rname)//&
877	" thres_n2ice = ",thres_n2ice
878
879	if (is_master)write(,)trim(rname)//&
880	"Diff of N2 albedo with thickness"
881	deltab=0. ! default value
882	call getin_p("deltab",deltab)
883	if (is_master)write(,)trim(rname)//&
884	" deltab = ",deltab
885
886	if (is_master)write(,)trim(rname)//&
887	"albedo N2 beta "
888	alb_n2b=0.7 ! default value
889	call getin_p("alb_n2b",alb_n2b)
890	if (is_master)write(,)trim(rname)//&
891	" alb_n2b = ",alb_n2b
892
893	if (is_master)write(,)trim(rname)//&
894	"albedo N2 alpha "
895	alb_n2a=0.7 ! default value
896	call getin_p("alb_n2a",alb_n2a)
897	if (is_master)write(,)trim(rname)//&
898	" alb_n2a = ",alb_n2a
899
900	if (is_master)write(,)trim(rname)//&
901	"emis N2 beta "
902	emis_n2b=0.7 ! default value
903	call getin_p("emis_n2b",emis_n2b)
904	if (is_master)write(,)trim(rname)//&
905	" emis_n2b = ",emis_n2b
906
907	if (is_master)write(,)trim(rname)//&
908	"emis N2 alpha "
909	emis_n2a=0.7 ! default value
910	call getin_p("emis_n2a",emis_n2a)
911	if (is_master)write(,)trim(rname)//&
912	" emis_n2a = ",emis_n2a
913
914	!********* CH4 *******************************
915
916	if (is_master)write(,)trim(rname)//&
917	"Mode for changing CH4 albedo"
918	mode_ch4=0 ! default value
919	call getin_p("mode_ch4",mode_ch4)
920	if (is_master)write(,)trim(rname)//&
921	" mode_ch4 = ",mode_ch4
922	feedback_met=0 ! default value
923	call getin_p("feedback_met",feedback_met)
924	if (is_master)write(,)trim(rname)//&
925	" feedback_met = ",feedback_met
926	thres_ch4ice=1. ! default value, mm
927	call getin_p("thres_ch4ice",thres_ch4ice)
928	if (is_master)write(,)trim(rname)//&
929	" thres_ch4ice = ",thres_ch4ice
930	fdch4_latn=-1.
931	fdch4_lats=0.
932	fdch4_lone=0.
933	fdch4_lonw=-1.
934	fdch4_depalb=0.5
935	fdch4_finalb=0.95
936	fdch4_maxalb=0.99
937	fdch4_ampl=200.
938	fdch4_maxice=100.
939	call getin_p("fdch4_latn",fdch4_latn)
940	call getin_p("fdch4_lats",fdch4_lats)
941	call getin_p("fdch4_lone",fdch4_lone)
942	call getin_p("fdch4_lonw",fdch4_lonw)
943	call getin_p("fdch4_depalb",fdch4_depalb)
944	call getin_p("fdch4_finalb",fdch4_finalb)
945	call getin_p("fdch4_maxalb",fdch4_maxalb)
946	call getin_p("fdch4_maxice",fdch4_maxice)
947	call getin_p("fdch4_ampl",fdch4_ampl)
948	if (is_master)write(,)trim(rname)//&
949	" Values for albedo feedback = ",fdch4_latn,&
950	fdch4_lats,fdch4_lone,fdch4_lonw,fdch4_depalb,&
951	fdch4_finalb,fdch4_maxalb,fdch4_maxice,fdch4_ampl
952
953	if (is_master)write(,)trim(rname)//&
954	"Latitude for diff albedo methane"
955	metlateq=25. ! default value
956	call getin_p("metlateq",metlateq)
957	if (is_master)write(,)trim(rname)//&
958	" metlateq = ",metlateq
959
960	if (is_master)write(,)trim(rname)//&
961	"Ls1 and Ls2 for change of ch4 albedo"
962	metls1=-1. ! default value
963	metls2=-2. ! default value
964	call getin_p("metls1",metls1)
965	call getin_p("metls2",metls2)
966
967	if (is_master)write(,)trim(rname)//&
968	"albedo CH4 "
969	alb_ch4=0.5 ! default value
970	call getin_p("alb_ch4",alb_ch4)
971	if (is_master)write(,)trim(rname)//&
972	" alb_ch4 = ",alb_ch4
973
974	if (is_master)write(,)trim(rname)//&
975	"albedo equatorial CH4 "
976	alb_ch4_eq=alb_ch4 ! default value
977	call getin_p("alb_ch4_eq",alb_ch4_eq)
978	if (is_master)write(,)trim(rname)//&
979	" alb_ch4_eq = ",alb_ch4_eq
980
981	if (is_master)write(,)trim(rname)//&
982	"albedo south hemis CH4 "
983	alb_ch4_s=alb_ch4 ! default value
984	call getin_p("alb_ch4_s",alb_ch4_s)
985	if (is_master)write(,)trim(rname)//&
986	" alb_ch4_s = ",alb_ch4_s
987
988	if (is_master)write(,)trim(rname)//&
989	"emis CH4 "
990	emis_ch4=0.5 ! default value
991	call getin_p("emis_ch4",emis_ch4)
992	if (is_master)write(,)trim(rname)//&
993	" emis_ch4 = ",emis_ch4
994
995	if (is_master)write(,)trim(rname)//&
996	"CH4 lag for n2 sublimation limitation"
997	ch4lag=.false. ! default value
998	latlag=-90. ! default value
999	vmrlag=1. ! default value
1000	call getin_p("ch4lag",ch4lag)
1001	call getin_p("latlag",latlag)
1002	call getin_p("vmrlag",vmrlag)
1003	if (is_master)write(,)trim(rname)//&
1004	" ch4lag = ",ch4lag
1005	if (is_master)write(,)trim(rname)//&
1006	" latlag = ",latlag
1007	if (is_master)write(,)trim(rname)//&
1008	" vmrlag = ",vmrlag
1009
1010	if (is_master)write(,)trim(rname)//&
1011	"Max temperature for surface ?"
1012	tsurfmax=.false. ! default value
1013	albmin_ch4=0.3 ! default value
1014	call getin_p("tsurfmax",tsurfmax)
1015	call getin_p("albmin_ch4",albmin_ch4)
1016	if (is_master)write(,)trim(rname)//&
1017	" tsurfmax = ",tsurfmax
1018	if (is_master)write(,)trim(rname)//&
1019	" albmin_ch4 = ",albmin_ch4
1020
1021
1022	if (is_master)write(,)trim(rname)//&
1023	"fixed gaseous CH4 mixing ratio for rad transfer ?"
1024	ch4fix=.false. ! default value
1025	call getin_p("ch4fix",ch4fix)
1026	if (is_master)write(,)trim(rname)//&
1027	" ch4fix = ",ch4fix
1028	if (is_master)write(,)trim(rname)//&
1029	"fixed gaseous CH4 mixing ratio for rad transfer ?"
1030	vmrch4fix=0.5 ! default value
1031	call getin_p("vmrch4fix",vmrch4fix)
1032	if (is_master)write(,)trim(rname)//&
1033	" vmrch4fix = ",vmrch4fix
1034	vmrch4_proffix=.false. ! default value
1035	call getin_p("vmrch4_proffix",vmrch4_proffix)
1036	if (is_master)write(,)trim(rname)//&
1037	" vmrch4_proffix = ",vmrch4_proffix
1038
1039
1040	!********* CO *******************************
1041
1042	if (is_master)write(,)trim(rname)//&
1043	"albedo CO "
1044	alb_co=0.4 ! default value
1045	call getin_p("alb_co",alb_co)
1046	if (is_master)write(,)trim(rname)//&
1047	" alb_co = ",alb_co
1048	thres_coice=1. ! default value, mm
1049	call getin_p("thres_coice",thres_coice)
1050	if (is_master)write(,)trim(rname)//&
1051	" thres_coice = ",thres_coice
1052
1053	if (is_master)write(,)trim(rname)//&
1054	"emis CO "
1055	emis_co=0.4 ! default value
1056	call getin_p("emis_co",emis_co)
1057	if (is_master)write(,)trim(rname)//&
1058	" emis_co = ",emis_co
1059
1060	!********* THOLINS *******************************
1061	if (is_master)write(,)trim(rname)//&
1062	"Mode for changing tholins albedo/emis"
1063	mode_tholins=0 ! default value
1064	call getin_p("mode_tholins",mode_tholins)
1065	if (is_master)write(,)trim(rname)//&
1066	" mode_tholins = ",mode_tholins
1067
1068	if (is_master)write(,)trim(rname)//&
1069	"albedo tho "
1070	alb_tho=0.1 ! default value
1071	call getin_p("alb_tho",alb_tho)
1072	if (is_master)write(,)trim(rname)//&
1073	" alb_tho = ",alb_tho
1074
1075	if (is_master)write(,)trim(rname)//&
1076	"albedo tho eq"
1077	alb_tho_eq=0.1 ! default value
1078	call getin_p("alb_tho_eq",alb_tho_eq)
1079	if (is_master)write(,)trim(rname)//&
1080	" alb_tho_eq = ",alb_tho_eq
1081
1082	if (is_master)write(,)trim(rname)//&
1083	"emis tholins "
1084	emis_tho=1. ! default value
1085	call getin_p("emis_tho",emis_tho)
1086	if (is_master)write(,)trim(rname)//&
1087	" emis_tho = ",emis_tho
1088
1089	if (is_master)write(,)trim(rname)//&
1090	"emis tholins eq"
1091	emis_tho_eq=1. ! default value
1092	call getin_p("emis_tho_eq",emis_tho_eq)
1093	if (is_master)write(,)trim(rname)//&
1094	" emis_tho_eq = ",emis_tho_eq
1095
1096	if (is_master)write(,)trim(rname)//&
1097	"Latitude for diff albedo tholins"
1098	tholateq=25. ! default value
1099	call getin_p("tholateq",tholateq)
1100	if (is_master)write(,)trim(rname)//&
1101	" tholateq = ",tholateq
1102	tholatn=-1.
1103	tholats=0.
1104	tholone=0.
1105	tholonw=-1.
1106	alb_tho_spe=0.1 ! default value
1107	emis_tho_spe=1. ! default value
1108	call getin_p(" tholatn",tholatn)
1109	call getin_p(" tholats",tholats)
1110	call getin_p(" tholone",tholone)
1111	call getin_p(" tholonw",tholonw)
1112	if (is_master)write(,)trim(rname)//&
1113	" Values for special tholins albedo = ",tholatn,&
1114	tholats,tholone,tholonw,alb_tho_spe,emis_tho_spe
1115
1116	if (is_master)write(,)trim(rname)//&
1117	"Specific albedo"
1118	spelon1=-180. ! default value
1119	spelon2=180. ! default value
1120	spelat1=-90. ! default value
1121	spelat2=90. ! default value
1122	specalb=.false. ! default value
1123	if (is_master)write(,)trim(rname)//&
1124	"albedo/emis spe "
1125	albspe=0.1 ! default value
1126	emispe=1. ! default value
1127	call getin_p("spelon1",spelon1)
1128	call getin_p("spelon2",spelon2)
1129	call getin_p("spelat1",spelat1)
1130	call getin_p("spelat2",spelat2)
1131	call getin_p("specalb",specalb)
1132	call getin_p("albspe",albspe)
1133	call getin_p("emispe",emispe)
1134
1135	if (is_master)write(,)trim(rname)//&
1136	" specific = ",specalb
1137
1138	!******************** TI ***************************
1139
1140	if (is_master)write(,)trim(rname)//&
1141	"Change TI with time"
1142	changeti=.false. ! default value
1143	call getin_p("changeti",changeti)
1144	if (is_master)write(,)trim(rname)//&
1145	" changeti = ",changeti
1146	changetid=.false. ! default value for diurnal TI
1147	call getin_p("changetid",changetid)
1148	if (is_master)write(,)trim(rname)//&
1149	" changetid = ",changetid
1150
1151	if (is_master)write(,)trim(rname)//&
1152	"IT N2 "
1153	ITN2=800. ! default value
1154	call getin_p("ITN2",ITN2)
1155	if (is_master)write(,)trim(rname)//&
1156	" ITN2 = ",ITN2
1157	ITN2d=20. ! default value
1158	call getin_p("ITN2d",ITN2d)
1159	if (is_master)write(,)trim(rname)//&
1160	" ITN2d = ",ITN2d
1161
1162	if (is_master)write(,)trim(rname)//&
1163	"IT CH4"
1164	ITCH4=800. ! default value
1165	call getin_p("ITCH4",ITCH4)
1166	if (is_master)write(,)trim(rname)//&
1167	" ITCH4 = ",ITCH4
1168	ITCH4d=20. ! default value
1169	call getin_p("ITCH4d",ITCH4d)
1170	if (is_master)write(,)trim(rname)//&
1171	" ITCH4d = ",ITCH4d
1172
1173	if (is_master)write(,)trim(rname)//&
1174	"IT H2O"
1175	ITH2O=800. ! default value
1176	call getin_p("ITH2O",ITH2O)
1177	if (is_master)write(,)trim(rname)//&
1178	" ITH2O = ",ITH2O
1179	ITH2Od=20. ! default value
1180	call getin_p("ITH2Od",ITH2Od)
1181	if (is_master)write(,)trim(rname)//&
1182	" ITH2Od = ",ITH2Od
1183
1184	!************ COOLING *************
1185
1186	alpha_top=5e-11 ! default value
1187	call getin_p("alpha_top",alpha_top)
1188	if (is_master)write(,)trim(rname)//&
1189	" alpha_top = ",alpha_top
1190	pref=0.02 ! default value
1191	call getin_p("pref",pref)
1192	if (is_master)write(,)trim(rname)//&
1193	" pref = ",pref
1194	deltap=0.1 ! default value
1195	call getin_p("deltap",deltap)
1196	if (is_master)write(,)trim(rname)//&
1197	" deltap = ",deltap
1198
1199	!=================================
1200	! TWOLAY scheme and NH3 cloudare left for retrocompatibility only,
1201	! You should now use N-LAYER scheme (see below).
1202
1203	if (is_master) write(,)trim(rname)//&
1204	": Radiatively active two-layer aerosols?"
1205	aeroback2lay=.false. ! default value
1206	call getin_p("aeroback2lay",aeroback2lay)
1207	if (is_master) write(,)trim(rname)//": aeroback2lay = ",aeroback2lay
1208
1209	if (aeroback2lay.and.is_master) then
1210	print*,'Warning : The TWOLAY AEROSOL scheme is deprecated and buggy...'
1211	print*,'You should use the generic n-layer scheme (see aeronlay).'
1212	endif
1213
1214	if (is_master) write(,)trim(rname)//": Radiatively active ammonia cloud?"
1215	aeronh3=.false. ! default value
1216	call getin_p("aeronh3",aeronh3)
1217	if (is_master) write(,)trim(rname)//": aeronh3 = ",aeronh3
1218
1219	if (aeronh3.and.is_master) then
1220	print*,'Warning : You are using specific NH3 cloud scheme ...'
1221	print*,'You should use the generic n-layer scheme (see aeronlay).'
1222	endif
1223
1224	if (is_master) write(,)trim(rname)//&
1225	": TWOLAY AEROSOL: total optical depth "//&
1226	"in the tropospheric layer (visible)"
1227	obs_tau_col_tropo=8.D0
1228	call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
1229	if (is_master) write(,)trim(rname)//&
1230	": obs_tau_col_tropo = ",obs_tau_col_tropo
1231
1232	if (is_master) write(,)trim(rname)//&
1233	": TWOLAY AEROSOL: total optical depth "//&
1234	"in the stratospheric layer (visible)"
1235	obs_tau_col_strato=0.08D0
1236	call getin_p("obs_tau_col_strato",obs_tau_col_strato)
1237	if (is_master) write(,)trim(rname)//&
1238	": obs_tau_col_strato = ",obs_tau_col_strato
1239
1240	if (is_master) write(,)trim(rname)//&
1241	": TWOLAY AEROSOL: optprop_back2lay_vis?"
1242	optprop_back2lay_vis = 'optprop_saturn_vis_n20.dat'
1243	call getin_p("optprop_back2lay_vis",optprop_back2lay_vis)
1244	if (is_master) write(,)trim(rname)//&
1245	": optprop_back2lay_vis = ",trim(optprop_back2lay_vis)
1246
1247	if (is_master) write(,)trim(rname)//&
1248	": TWOLAY AEROSOL: optprop_back2lay_ir?"
1249	optprop_back2lay_ir = 'optprop_saturn_ir_n20.dat'
1250	call getin_p("optprop_back2lay_ir",optprop_back2lay_ir)
1251	if (is_master) write(,)trim(rname)//&
1252	": optprop_back2lay_ir = ",trim(optprop_back2lay_ir)
1253
1254	if (is_master) write(,)trim(rname)//&
1255	": TWOLAY AEROSOL: pres_bottom_tropo? in pa"
1256	pres_bottom_tropo=66000.0
1257	call getin_p("pres_bottom_tropo",pres_bottom_tropo)
1258	if (is_master) write(,)trim(rname)//&
1259	": pres_bottom_tropo = ",pres_bottom_tropo
1260
1261	if (is_master) write(,)trim(rname)//&
1262	": TWOLAY AEROSOL: pres_top_tropo? in pa"
1263	pres_top_tropo=18000.0
1264	call getin_p("pres_top_tropo",pres_top_tropo)
1265	if (is_master) write(,)trim(rname)//&
1266	": pres_top_tropo = ",pres_top_tropo
1267
1268	if (is_master) write(,)trim(rname)//&
1269	": TWOLAY AEROSOL: pres_bottom_strato? in pa"
1270	pres_bottom_strato=2000.0
1271	call getin_p("pres_bottom_strato",pres_bottom_strato)
1272	if (is_master) write(,)trim(rname)//&
1273	": pres_bottom_strato = ",pres_bottom_strato
1274
1275	! Sanity check
1276	if (pres_bottom_strato .gt. pres_top_tropo) then
1277	if(is_master) then
1278	print*,'Error : TWOLAY AEROSOL, Please ensure that in callphys.def'
1279	print*,'you have pres_top_tropo > pres_bottom_strato !'
1280	endif
1281	call abort_physic(rname," pres_top_tropo > pres_bottom_strato!",1)
1282	endif
1283
1284	if (is_master) write(,)trim(rname)//&
1285	": TWOLAY AEROSOL: pres_top_strato? in pa"
1286	pres_top_strato=100.0
1287	call getin_p("pres_top_strato",pres_top_strato)
1288	if (is_master) write(,)trim(rname)//&
1289	": pres_top_strato = ",pres_top_strato
1290
1291	if (is_master) write(,)trim(rname)//&
1292	": TWOLAY AEROSOL: particle size in the ", &
1293	"tropospheric layer, in meters"
1294	size_tropo=2.e-6
1295	call getin_p("size_tropo",size_tropo)
1296	if (is_master) write(,)trim(rname)//": size_tropo = ",size_tropo
1297
1298	if (is_master) write(,)trim(rname)//&
1299	": TWOLAY AEROSOL: particle size in the ", &
1300	"stratospheric layer, in meters"
1301	size_strato=1.e-7
1302	call getin_p("size_strato",size_strato)
1303	if (is_master) write(,)trim(rname)//": size_strato = ",size_strato
1304
1305	if (is_master) write(,)trim(rname)//&
1306	": NH3 (thin) cloud: total optical depth"
1307	tau_nh3_cloud=7.D0
1308	call getin_p("tau_nh3_cloud",tau_nh3_cloud)
1309	if (is_master) write(,)trim(rname)//": tau_nh3_cloud = ",tau_nh3_cloud
1310
1311	if (is_master) write(,)trim(rname)//": NH3 (thin) cloud pressure level"
1312	pres_nh3_cloud=7.D0
1313	call getin_p("pres_nh3_cloud",pres_nh3_cloud)
1314	if (is_master) write(,)trim(rname)//": pres_nh3_cloud = ",pres_nh3_cloud
1315
1316	if (is_master) write(,)trim(rname)//": NH3 (thin) cloud: particle sizes"
1317	size_nh3_cloud=3.e-6
1318	call getin_p("size_nh3_cloud",size_nh3_cloud)
1319	if (is_master) write(,)trim(rname)//": size_nh3_cloud = ",size_nh3_cloud
1320
1321	!=================================
1322	! Generic N-LAYER aerosol scheme
1323
1324	if (is_master) write(,)trim(rname)//&
1325	": Radiatively active generic n-layer aerosols?"
1326	aeronlay=.false. ! default value
1327	call getin_p("aeronlay",aeronlay)
1328	if (is_master) write(,)trim(rname)//": aeronlay = ",aeronlay
1329
1330	if (is_master) write(,)trim(rname)//&
1331	": Number of generic aerosols layers?"
1332	nlayaero=1 ! default value
1333	call getin_p("nlayaero",nlayaero)
1334	! Avoid to allocate arrays of size 0
1335	if (aeronlay .and. nlayaero.lt.1) then
1336	if (is_master) then
1337	print*, " You are trying to set no generic aerosols..."
1338	print*, " Set aeronlay=.false. instead ! I abort."
1339	endif
1340	call abort_physic(rname,"no generic aerosols...",1)
1341	endif
1342	if (.not. aeronlay) nlayaero=1
1343	if (is_master) write(,)trim(rname)//": nlayaero = ",nlayaero
1344
1345	! This is necessary, we just set the number of aerosol layers
1346	IF(.NOT.ALLOCATED(aeronlay_tauref)) ALLOCATE(aeronlay_tauref(nlayaero))
1347	IF(.NOT.ALLOCATED(aeronlay_lamref)) ALLOCATE(aeronlay_lamref(nlayaero))
1348	IF(.NOT.ALLOCATED(aeronlay_choice)) ALLOCATE(aeronlay_choice(nlayaero))
1349	IF(.NOT.ALLOCATED(aeronlay_pbot)) ALLOCATE(aeronlay_pbot(nlayaero))
1350	IF(.NOT.ALLOCATED(aeronlay_ptop)) ALLOCATE(aeronlay_ptop(nlayaero))
1351	IF(.NOT.ALLOCATED(aeronlay_sclhght)) ALLOCATE(aeronlay_sclhght(nlayaero))
1352	IF(.NOT.ALLOCATED(aeronlay_size)) ALLOCATE(aeronlay_size(nlayaero))
1353	IF(.NOT.ALLOCATED(aeronlay_nueff)) ALLOCATE(aeronlay_nueff(nlayaero))
1354	IF(.NOT.ALLOCATED(optprop_aeronlay_ir)) ALLOCATE(optprop_aeronlay_ir(nlayaero))
1355	IF(.NOT.ALLOCATED(optprop_aeronlay_vis)) ALLOCATE(optprop_aeronlay_vis(nlayaero))
1356
1357	if (is_master) write(,)trim(rname)//&
1358	": Generic n-layer aerosols: Optical depth at reference wavelenght"
1359	aeronlay_tauref=1.0E-1
1360	call getin_p("aeronlay_tauref",aeronlay_tauref)
1361	if (is_master) write(,)trim(rname)//&
1362	": aeronlay_tauref = ",aeronlay_tauref
1363
1364	if (is_master) write(,)trim(rname)//&
1365	": Generic n-layer aerosols: Reference wavelenght for optical depths (m)"
1366	aeronlay_lamref=0.6E-6
1367	call getin_p("aeronlay_lamref",aeronlay_lamref)
1368	if (is_master) write(,)trim(rname)//&
1369	": aeronlay_lamref = ",aeronlay_lamref
1370
1371	if (is_master) then
1372	write(,)trim(rname)//&
1373	": Generic n-layer aerosols: Vertical profile choice : "
1374	write(,)trim(rname)//&
1375	" (1) Tau btwn ptop and pbot follows atm. scale height"
1376	write(,)trim(rname)//&
1377	" or (2) Tau above pbot follows its own scale height"
1378	endif
1379	aeronlay_choice=1
1380	call getin_p("aeronlay_choice",aeronlay_choice)
1381	if (is_master) write(,)trim(rname)//&
1382	": aeronlay_choice = ",aeronlay_choice
1383
1384	if (is_master) write(,)trim(rname)//&
1385	": Generic n-layer aerosols: bottom pressures (Pa)"
1386	aeronlay_pbot=2000.0
1387	call getin_p("aeronlay_pbot",aeronlay_pbot)
1388	if (is_master) write(,)trim(rname)//": aeronlay_pbot = ",aeronlay_pbot
1389
1390	if (is_master) write(,)trim(rname)//&
1391	": Generic n-layer aerosols: (if choice=1) Top pressures (Pa) "
1392	aeronlay_ptop=300000.0
1393	call getin_p("aeronlay_ptop",aeronlay_ptop)
1394	if (is_master) write(,)trim(rname)//": aeronlay_ptop = ",aeronlay_ptop
1395
1396	if (is_master) write(,)trim(rname)//&
1397	": Generic n-layer aerosols: (if choice=2) Scale height / atm. scale height"
1398	aeronlay_sclhght=0.2
1399	call getin_p("aeronlay_sclhght",aeronlay_sclhght)
1400	if (is_master) write(,)trim(rname)//&
1401	": aeronlay_sclhght = ",aeronlay_sclhght
1402
1403	if (is_master) write(,)trim(rname)//&
1404	": Generic n-layer aerosols: particles effective radii (m)"
1405	aeronlay_size=1.e-6
1406	call getin_p("aeronlay_size",aeronlay_size)
1407	if (is_master) write(,)trim(rname)//": aeronlay_size = ",aeronlay_size
1408
1409	if (is_master) write(,)trim(rname)//&
1410	": Generic n-layer aerosols: particles radii effective variance"
1411	aeronlay_nueff=0.1
1412	call getin_p("aeronlay_nueff",aeronlay_nueff)
1413	if (is_master) write(,)trim(rname)//": aeronlay_nueff = ",aeronlay_nueff
1414
1415	if (is_master) write(,)trim(rname)//&
1416	": Generic n-layer aerosols: VIS optical properties file"
1417	optprop_aeronlay_vis = 'optprop_saturn_vis_n20.dat'
1418	call getin_p("optprop_aeronlay_vis",optprop_aeronlay_vis)
1419	if (is_master) write(,)trim(rname)//&
1420	": optprop_aeronlay_vis = ",optprop_aeronlay_vis
1421
1422	if (is_master) write(,)trim(rname)//&
1423	": Generic n-layer aerosols: IR optical properties file"
1424	optprop_aeronlay_ir = 'optprop_saturn_ir_n20.dat'
1425	call getin_p("optprop_aeronlay_ir",optprop_aeronlay_ir)
1426	if (is_master) write(,)trim(rname)//&
1427	": optprop_aeronlay_ir = ",optprop_aeronlay_ir
1428
1429
1430	!=================================
1431
1432	if (is_master) write(,)trim(rname)//&
1433	": Is the variable gas species radiatively active?"
1434	Tstrat=167.0
1435	varactive=.false.
1436	call getin_p("varactive",varactive)
1437	if (is_master) write(,)trim(rname)//": varactive = ",varactive
1438
1439	if (is_master) write(,)trim(rname)//&
1440	": Is the variable gas species distribution set?"
1441	varfixed=.false.
1442	call getin_p("varfixed",varfixed)
1443	if (is_master) write(,)trim(rname)//": varfixed = ",varfixed
1444
1445	if (is_master) write(,)trim(rname)//&
1446	": What is the saturation % of the variable species?"
1447	satval=0.8
1448	call getin_p("satval",satval)
1449	if (is_master) write(,)trim(rname)//": satval = ",satval
1450
1451
1452	! Test of incompatibility:
1453	! if varactive, then varfixed should be false
1454	if (varactive.and.varfixed) then
1455	call abort_physic(rname,'if varactive, varfixed must be OFF!',1)
1456	endif
1457
1458	if (is_master) write(,)trim(rname)//": Gravitationnal sedimentation ?"
1459	sedimentation=.false. ! default value
1460	call getin_p("sedimentation",sedimentation)
1461	if (is_master) write(,)trim(rname)//": sedimentation = ",sedimentation
1462
1463	if (is_master) write(,)trim(rname)//": Generic Condensation of tracers ?"
1464	generic_condensation=.false. !default value
1465	call getin_p("generic_condensation",generic_condensation)
1466	if (is_master) write(,)trim(rname)//": generic_condensation = ",generic_condensation
1467
1468	if (is_master) write(,)trim(rname)//": Spectral Dependant albedo ?"
1469	albedo_spectral_mode=.false. ! default value
1470	call getin_p("albedo_spectral_mode",albedo_spectral_mode)
1471	if (is_master) write(,)trim(rname)//&
1472	": albedo_spectral_mode = ",albedo_spectral_mode
1473
1474	if (is_master) then
1475	write(,)trim(rname)//": Snow albedo ?"
1476	write(,)trim(rname)//": If albedo_spectral_mode=.true., then this "
1477	write(,)trim(rname)//": corresponds to the 0.5 microns snow albedo."
1478	endif
1479	albedosnow=0.5 ! default value
1480	call getin_p("albedosnow",albedosnow)
1481	if (is_master) write(,)trim(rname)//": albedosnow = ",albedosnow
1482
1483	if (is_master) write(,)trim(rname)//": N2 ice albedo ?"
1484	albedon2ice=0.5 ! default value
1485	call getin_p("albedon2ice",albedon2ice)
1486	if (is_master) write(,)trim(rname)//": albedon2ice = ",albedon2ice
1487
1488	if (is_master) write(,)trim(rname)//": Maximum ice thickness ?"
1489	maxicethick=2.0 ! default value
1490	call getin_p("maxicethick",maxicethick)
1491	if (is_master) write(,)trim(rname)//": maxicethick = ",maxicethick
1492
1493	if (is_master) write(,)trim(rname)//": Minimum eddy mix coeff in 1D ?"
1494	kmixmin=1.0e-2 ! default value
1495	call getin_p("kmixmin",kmixmin)
1496	if (is_master) write(,)trim(rname)//": kmixmin = ",kmixmin
1497
1498	if (is_master) write(,)'Predefined Cp from dynamics is ',cpp,'J kg^-1 K^-1'
1499	if (is_master) write(,)'Predefined Mg from dynamics is ',mugaz,'amu'
1500
1501	force_cpp=.false. ! default value
1502	call getin_p("force_cpp",force_cpp)
1503	if (force_cpp) then
1504	if (is_master) write(,)trim(rname)//": force_cpp = ",force_cpp
1505	if (is_master) write(,)trim(rname)//": force_cpp is deprecated.",&
1506	"Set cpp_mugaz_mode=1 in callfis to emulate force_cpp=.true."
1507	call abort_physic(rname,"Anyway, you need to set force_cpp=.false. to continue.",1)
1508	endif
1509
1510	if (is_master) write(,)trim(rname)//&
1511	": where do you want your cpp/mugaz value to come from?",&
1512	"=> 0: dynamics (3d), 1: forced in callfis (1d), 2: computed from gases.def (1d)?"
1513	cpp_mugaz_mode = 0 ! default value
1514	call getin_p("cpp_mugaz_mode",cpp_mugaz_mode)
1515	if (is_master) write(,)trim(rname)//": cpp_mugaz_mode = ",cpp_mugaz_mode
1516
1517
1518	! Test of incompatibility:
1519
1520	if ((.not.tracer).and.(haze)) then
1521	call abort_physic(rname, 'if haze are on, tracers must be on!', 1)
1522	endif
1523	if (haze_proffix.and.sedimentation) then
1524	call abort_physic(rname, 'if haze profile is set, sedimentation must be deactivated', 1)
1525	endif
1526	if (callmolvis.and..not.callconduct) then
1527	call abort_physic(rname, 'if callmolvis is set, callconduct must be true', 1)
1528	endif
1529	if (glaflow.and..not.fast) then
1530	call abort_physic(rname, 'if glaflow is set, fast must be true', 1)
1531	endif
1532	if (paleo.and..not.fast) then
1533	call abort_physic(rname, 'if paleo is set, fast must be true', 1)
1534	endif
1535	if ((haze_proffix.or.haze_radproffix).and..not.aerohaze) then
1536	call abort_physic(rname, 'for now, haze/rad proffix only works w aerohaze=T', 1)
1537	endif
1538	if (carbox.and.condcosurf.and.no_n2frost) then
1539	call abort_physic(rname, "CO surface condensation and no_n2frost are both active which may not be relevant", 1)
1540	end if
1541
1542	if ((cpp_mugaz_mode >= 1).and.(is_master).and.(ngrid>1)) then
1543	write(,)' !!!!! Be aware that having different values for cpp and mugaz in the dynamics and physics'
1544	write(,)' in 3D can result in very pathological behavior. You have been warned !!!!!'
1545	else if ((cpp_mugaz_mode == 0).and.(is_master).and.(ngrid==1)) then
1546	! for this specific 1D error we will remove run.def before aborting JL22
1547	call system("rm -rf run.def")
1548	call abort_physic(rname,"cpp_mugaz_mode must be >= 1 in 1d",1)
1549	endif
1550
1551	if (cpp_mugaz_mode == 1) then
1552	mugaz = -99999.
1553	if (is_master) write(,)trim(rname)//&
1554	": MEAN MOLECULAR MASS in g mol-1 ?"
1555	call getin_p("mugaz",mugaz)
1556	IF (mugaz.eq.-99999.) THEN
1557	call abort_physic(rname,"mugaz must be set if cpp_mugaz_mode = 1",1)
1558	ENDIF
1559	cpp = -99999.
1560	if (is_master) write(,)trim(rname)//&
1561	": SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
1562	call getin_p("cpp",cpp)
1563	IF (cpp.eq.-99999.) THEN
1564	PRINT *, "cpp must be set if cpp_mugaz_mode = 1"
1565	STOP
1566	ENDIF
1567	if (is_master) write(,)'New Cp from callfis is ',cpp,'J kg^-1 K^-1'
1568	if (is_master) write(,)'New Mg from callfis is ',mugaz,'amu'
1569
1570	endif ! of if (cpp_mugaz_mode == 1)
1571	call su_gases
1572	call calc_cpp_mugaz
1573
1574	if (is_master) then
1575	PRINT*,'--------------------------------------------'
1576	PRINT*
1577	PRINT*
1578	endif
1579	ELSE
1580	call abort_physic(rname,'Cannot read file callphys.def. Is it here ?',1)
1581	ENDIF ! of IF(iscallphys)
1582
1583	if (is_master) then
1584	PRINT*
1585	PRINT*,'inifis: daysec',daysec
1586	PRINT*
1587	PRINT*,'inifis: The radiative transfer is computed:'
1588	PRINT*,' each ',iradia,' physical time-step'
1589	PRINT,' or each ',iradiadtphys,' seconds'
1590	PRINT*
1591	endif
1592
1593	!-----------------------------------------------------------------------
1594	! Some more initialization:
1595	! ------------------------
1596
1597	! Initializations for comgeomfi_h
1598	#ifndef MESOSCALE
1599	totarea=SSUM(ngrid,parea,1)
1600	call planetwide_sumval(parea,totarea_planet)
1601
1602	!! those are defined in comdiurn_h.F90
1603	IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
1604	IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
1605	IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
1606	IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
1607
1608	DO ig=1,ngrid
1609	sinlat(ig)=sin(plat(ig))
1610	coslat(ig)=cos(plat(ig))
1611	sinlon(ig)=sin(plon(ig))
1612	coslon(ig)=cos(plon(ig))
1613	ENDDO
1614	#endif
1615	! initialize variables in radinc_h
1616	call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
1617
1618	! initialize variables and allocate arrays in radcommon_h
1619	call ini_radcommon_h(naerkind)
1620
1621	! allocate "comsoil_h" arrays
1622	call ini_comsoil_h(ngrid)
1623
1624	END SUBROUTINE inifis
1625
1626	END MODULE inifis_mod

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

Download in other formats:

Original Format