Context Navigation

inifis_mod.F90 @ 3474

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

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

Download in other formats:

Original Format