Context Navigation

inifis_mod.F90 @ 3542

Last change on this file since 3542 was 3539, checked in by tbertrand, 13 months ago
LMDZ.PLUTO: Update and Validation of the Volatile Transport Model (VTM, or "nogcm") Merging missing elements with Pluto.old TB
File size: 57.3 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 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	if (is_master)write(,)trim(rname)//&
672	"Production of haze ?"
673	haze=.false. ! default value
674	call getin_p("haze",haze)
675	if (is_master)write(,)trim(rname)//&
676	" haze = ",haze
677	hazeconservch4=.false. ! conservch4, by default value ch4 is photolyzed
678	call getin_p("hazeconservch4",hazeconservch4)
679	if (is_master)write(,)trim(rname)//&
680	"hazconservch4 = ",hazeconservch4
681	if (is_master)write(,)trim(rname)//&
682	"Production of haze (fast version) ?"
683	fasthaze=.false. ! default value
684	call getin_p("fasthaze",fasthaze)
685	if (is_master)write(,)trim(rname)//&
686	"fasthaze = ",fasthaze
687
688	if (is_master)write(,)trim(rname)//&
689	"Add source of haze ?"
690	source_haze=.false. ! default value
691	call getin_p("source_haze",source_haze)
692	if (is_master)write(,)trim(rname)//&
693	" source_haze = ",source_haze
694	mode_hs=0 ! mode haze source default value
695	call getin_p("mode_hs",mode_hs)
696	if (is_master)write(,)trim(rname)//&
697	" mode_hs",mode_hs
698	kfix=1 ! default value
699	call getin_p("kfix",kfix)
700	if (is_master)write(,)trim(rname)//&
701	"levels kfix",kfix
702	fracsource=0.1 ! default value
703	call getin_p("fracsource",fracsource)
704	if (is_master)write(,)trim(rname)//&
705	" fracsource",fracsource
706	latsource=30. ! default value
707	call getin_p("latsource",latsource)
708	if (is_master)write(,)trim(rname)//&
709	" latsource",latsource
710	lonsource=180. ! default value
711	call getin_p("lonsource",lonsource)
712	if (is_master)write(,)trim(rname)//&
713	" lonsource",lonsource
714
715	if (is_master)write(,)trim(rname)//&
716	"Radiatively active haze ?"
717	aerohaze=.false. ! default value
718	call getin_p("aerohaze",aerohaze)
719	if (is_master)write(,)trim(rname)//&
720	"aerohaze = ",aerohaze
721
722	if (is_master)write(,)trim(rname)//&
723	"Haze monomer radius ?"
724	rad_haze=10.e-9 ! default value
725	call getin_p("rad_haze",rad_haze)
726	if (is_master)write(,)trim(rname)//&
727	"rad_haze = ",rad_haze
728
729	if (is_master)write(,)trim(rname)//&
730	"fractal particle ?"
731	fractal=.false. ! default value
732	call getin_p("fractal",fractal)
733	if (is_master)write(,)trim(rname)//&
734	"fractal = ",fractal
735	nb_monomer=10 ! default value
736	call getin_p("nb_monomer",nb_monomer)
737	if (is_master)write(,)trim(rname)//&
738	"nb_monomer = ",nb_monomer
739
740	if (is_master)write(,)trim(rname)//&
741	"fixed gaseous CH4 mixing ratio for rad transfer ?"
742	ch4fix=.false. ! default value
743	call getin_p("ch4fix",ch4fix)
744	if (is_master)write(,)trim(rname)//&
745	" ch4fix = ",ch4fix
746	if (is_master)write(,)trim(rname)//&
747	"fixed gaseous CH4 mixing ratio for rad transfer ?"
748	vmrch4fix=0.5 ! default value
749	call getin_p("vmrch4fix",vmrch4fix)
750	if (is_master)write(,)trim(rname)//&
751	" vmrch4fix = ",vmrch4fix
752	vmrch4_proffix=.false. ! default value
753	call getin_p("vmrch4_proffix",vmrch4_proffix)
754	if (is_master)write(,)trim(rname)//&
755	" vmrch4_proffix = ",vmrch4_proffix
756
757	if (is_master)write(,)trim(rname)//&
758	"call specific cooling for radiative transfer ?"
759	cooling=.false. ! default value
760	call getin_p("cooling",cooling)
761	if (is_master)write(,)trim(rname)//&
762	" cooling = ",cooling
763
764	if (is_master)write(,)trim(rname)//&
765	"NLTE correction for methane heating rates?"
766	nlte=.false. ! default value
767	call getin_p("nlte",nlte)
768	if (is_master)write(,)trim(rname)//&
769	" nlte = ",nlte
770	strobel=.false. ! default value
771	call getin_p("strobel",strobel)
772	if (is_master)write(,)trim(rname)//&
773	" strobel = ",strobel
774
775	if (is_master)write(,)trim(rname)//&
776	"fixed radius profile from txt file ?"
777	haze_radproffix=.false. ! default value
778	call getin_p("haze_radproffix",haze_radproffix)
779	if (is_master)write(,)trim(rname)//&
780	"haze_radproffix = ",haze_radproffix
781	if (is_master)write(,)trim(rname)//&
782	"fixed MMR profile from txt file ?"
783	haze_proffix=.false. ! default value
784	call getin_p("haze_proffix",haze_proffix)
785	if (is_master)write(,)trim(rname)//&
786	"haze_proffix = ",haze_proffix
787
788	if (is_master)write(,)trim(rname)//&
789	"Number mix ratio of haze particles for co clouds:"
790	Nmix_co=100000. ! default value
791	call getin_p("Nmix_co",Nmix_co)
792	if (is_master)write(,)trim(rname)//&
793	" Nmix_co = ",Nmix_co
794
795	if (is_master)write(,)trim(rname)//&
796	"Number mix ratio of haze particles for ch4 clouds:"
797	Nmix_ch4=100000. ! default value
798	call getin_p("Nmix_ch4",Nmix_ch4)
799	if (is_master)write(,)trim(rname)//&
800	" Nmix_ch4 = ",Nmix_ch4
801
802	!***************************************************************
803	!! Surface properties
804
805	!********* N2 *******************************
806
807	if (is_master)write(,)trim(rname)//&
808	"Mode for changing N2 albedo"
809	mode_n2=0 ! default value
810	call getin_p("mode_n2",mode_n2)
811	if (is_master)write(,)trim(rname)//&
812	" mode_n2 = ",mode_n2
813	thres_n2ice=1. ! default value
814	call getin_p("thres_n2ice",thres_n2ice)
815	if (is_master)write(,)trim(rname)//&
816	" thres_n2ice = ",thres_n2ice
817
818	if (is_master)write(,)trim(rname)//&
819	"Diff of N2 albedo with thickness"
820	deltab=0. ! default value
821	call getin_p("deltab",deltab)
822	if (is_master)write(,)trim(rname)//&
823	" deltab = ",deltab
824
825	if (is_master)write(,)trim(rname)//&
826	"albedo N2 beta "
827	alb_n2b=0.7 ! default value
828	call getin_p("alb_n2b",alb_n2b)
829	if (is_master)write(,)trim(rname)//&
830	" alb_n2b = ",alb_n2b
831
832	if (is_master)write(,)trim(rname)//&
833	"albedo N2 alpha "
834	alb_n2a=0.7 ! default value
835	call getin_p("alb_n2a",alb_n2a)
836	if (is_master)write(,)trim(rname)//&
837	" alb_n2a = ",alb_n2a
838
839	if (is_master)write(,)trim(rname)//&
840	"emis N2 beta "
841	emis_n2b=0.7 ! default value
842	call getin_p("emis_n2b",emis_n2b)
843	if (is_master)write(,)trim(rname)//&
844	" emis_n2b = ",emis_n2b
845
846	if (is_master)write(,)trim(rname)//&
847	"emis N2 alpha "
848	emis_n2a=0.7 ! default value
849	call getin_p("emis_n2a",emis_n2a)
850	if (is_master)write(,)trim(rname)//&
851	" emis_n2a = ",emis_n2a
852
853	!********* CH4 *******************************
854
855	if (is_master)write(,)trim(rname)//&
856	"Mode for changing CH4 albedo"
857	mode_ch4=0 ! default value
858	call getin_p("mode_ch4",mode_ch4)
859	if (is_master)write(,)trim(rname)//&
860	" mode_ch4 = ",mode_ch4
861	feedback_met=0 ! default value
862	call getin_p("feedback_met",feedback_met)
863	if (is_master)write(,)trim(rname)//&
864	" feedback_met = ",feedback_met
865	thres_ch4ice=1. ! default value, mm
866	call getin_p("thres_ch4ice",thres_ch4ice)
867	if (is_master)write(,)trim(rname)//&
868	" thres_ch4ice = ",thres_ch4ice
869	fdch4_latn=-1.
870	fdch4_lats=0.
871	fdch4_lone=0.
872	fdch4_lonw=-1.
873	fdch4_depalb=0.5
874	fdch4_finalb=0.95
875	fdch4_maxalb=0.99
876	fdch4_ampl=200.
877	fdch4_maxice=100.
878	call getin_p("fdch4_latn",fdch4_latn)
879	call getin_p("fdch4_lats",fdch4_lats)
880	call getin_p("fdch4_lone",fdch4_lone)
881	call getin_p("fdch4_lonw",fdch4_lonw)
882	call getin_p("fdch4_depalb",fdch4_depalb)
883	call getin_p("fdch4_finalb",fdch4_finalb)
884	call getin_p("fdch4_maxalb",fdch4_maxalb)
885	call getin_p("fdch4_maxice",fdch4_maxice)
886	call getin_p("fdch4_ampl",fdch4_ampl)
887	if (is_master)write(,)trim(rname)//&
888	" Values for albedo feedback = ",fdch4_latn,&
889	fdch4_lats,fdch4_lone,fdch4_lonw,fdch4_depalb,&
890	fdch4_finalb,fdch4_maxalb,fdch4_maxice,fdch4_ampl
891
892	if (is_master)write(,)trim(rname)//&
893	"Latitude for diff albedo methane"
894	metlateq=25. ! default value
895	call getin_p("metlateq",metlateq)
896	if (is_master)write(,)trim(rname)//&
897	" metlateq = ",metlateq
898
899	if (is_master)write(,)trim(rname)//&
900	"Ls1 and Ls2 for change of ch4 albedo"
901	metls1=-1. ! default value
902	metls2=-2. ! default value
903	call getin_p("metls1",metls1)
904	call getin_p("metls2",metls2)
905
906	if (is_master)write(,)trim(rname)//&
907	"albedo CH4 "
908	alb_ch4=0.5 ! default value
909	call getin_p("alb_ch4",alb_ch4)
910	if (is_master)write(,)trim(rname)//&
911	" alb_ch4 = ",alb_ch4
912
913	if (is_master)write(,)trim(rname)//&
914	"albedo equatorial CH4 "
915	alb_ch4_eq=alb_ch4 ! default value
916	call getin_p("alb_ch4_eq",alb_ch4_eq)
917	if (is_master)write(,)trim(rname)//&
918	" alb_ch4_eq = ",alb_ch4_eq
919
920	if (is_master)write(,)trim(rname)//&
921	"albedo south hemis CH4 "
922	alb_ch4_s=alb_ch4 ! default value
923	call getin_p("alb_ch4_s",alb_ch4_s)
924	if (is_master)write(,)trim(rname)//&
925	" alb_ch4_s = ",alb_ch4_s
926
927	if (is_master)write(,)trim(rname)//&
928	"emis CH4 "
929	emis_ch4=0.5 ! default value
930	call getin_p("emis_ch4",emis_ch4)
931	if (is_master)write(,)trim(rname)//&
932	" emis_ch4 = ",emis_ch4
933
934	if (is_master)write(,)trim(rname)//&
935	"CH4 lag for n2 sublimation limitation"
936	ch4lag=.false. ! default value
937	latlag=-90. ! default value
938	vmrlag=1. ! default value
939	call getin_p("ch4lag",ch4lag)
940	call getin_p("latlag",latlag)
941	call getin_p("vmrlag",vmrlag)
942	if (is_master)write(,)trim(rname)//&
943	" ch4lag = ",ch4lag
944	if (is_master)write(,)trim(rname)//&
945	" latlag = ",latlag
946	if (is_master)write(,)trim(rname)//&
947	" vmrlag = ",vmrlag
948
949	if (is_master)write(,)trim(rname)//&
950	"Max temperature for surface ?"
951	tsurfmax=.false. ! default value
952	albmin_ch4=0.3 ! default value
953	call getin_p("tsurfmax",tsurfmax)
954	call getin_p("albmin_ch4",albmin_ch4)
955	if (is_master)write(,)trim(rname)//&
956	" tsurfmax = ",tsurfmax
957	if (is_master)write(,)trim(rname)//&
958	" albmin_ch4 = ",albmin_ch4
959
960
961	if (is_master)write(,)trim(rname)//&
962	"fixed gaseous CH4 mixing ratio for rad transfer ?"
963	ch4fix=.false. ! default value
964	call getin_p("ch4fix",ch4fix)
965	if (is_master)write(,)trim(rname)//&
966	" ch4fix = ",ch4fix
967	if (is_master)write(,)trim(rname)//&
968	"fixed gaseous CH4 mixing ratio for rad transfer ?"
969	vmrch4fix=0.5 ! default value
970	call getin_p("vmrch4fix",vmrch4fix)
971	if (is_master)write(,)trim(rname)//&
972	" vmrch4fix = ",vmrch4fix
973	vmrch4_proffix=.false. ! default value
974	call getin_p("vmrch4_proffix",vmrch4_proffix)
975	if (is_master)write(,)trim(rname)//&
976	" vmrch4_proffix = ",vmrch4_proffix
977
978
979	!********* CO *******************************
980
981	if (is_master)write(,)trim(rname)//&
982	"albedo CO "
983	alb_co=0.4 ! default value
984	call getin_p("alb_co",alb_co)
985	if (is_master)write(,)trim(rname)//&
986	" alb_co = ",alb_co
987	thres_coice=1. ! default value, mm
988	call getin_p("thres_coice",thres_coice)
989	if (is_master)write(,)trim(rname)//&
990	" thres_coice = ",thres_coice
991
992	if (is_master)write(,)trim(rname)//&
993	"emis CO "
994	emis_co=0.4 ! default value
995	call getin_p("emis_co",emis_co)
996	if (is_master)write(,)trim(rname)//&
997	" emis_co = ",emis_co
998
999	!********* THOLINS *******************************
1000	if (is_master)write(,)trim(rname)//&
1001	"Mode for changing tholins albedo/emis"
1002	mode_tholins=0 ! default value
1003	call getin_p("mode_tholins",mode_tholins)
1004	if (is_master)write(,)trim(rname)//&
1005	" mode_tholins = ",mode_tholins
1006
1007	if (is_master)write(,)trim(rname)//&
1008	"albedo tho "
1009	alb_tho=0.1 ! default value
1010	call getin_p("alb_tho",alb_tho)
1011	if (is_master)write(,)trim(rname)//&
1012	" alb_tho = ",alb_tho
1013
1014	if (is_master)write(,)trim(rname)//&
1015	"albedo tho eq"
1016	alb_tho_eq=0.1 ! default value
1017	call getin_p("alb_tho_eq",alb_tho_eq)
1018	if (is_master)write(,)trim(rname)//&
1019	" alb_tho_eq = ",alb_tho_eq
1020
1021	if (is_master)write(,)trim(rname)//&
1022	"emis tholins "
1023	emis_tho=1. ! default value
1024	call getin_p("emis_tho",emis_tho)
1025	if (is_master)write(,)trim(rname)//&
1026	" emis_tho = ",emis_tho
1027
1028	if (is_master)write(,)trim(rname)//&
1029	"emis tholins eq"
1030	emis_tho_eq=1. ! default value
1031	call getin_p("emis_tho_eq",emis_tho_eq)
1032	if (is_master)write(,)trim(rname)//&
1033	" emis_tho_eq = ",emis_tho_eq
1034
1035	if (is_master)write(,)trim(rname)//&
1036	"Latitude for diff albedo tholins"
1037	tholateq=25. ! default value
1038	call getin_p("tholateq",tholateq)
1039	if (is_master)write(,)trim(rname)//&
1040	" tholateq = ",tholateq
1041	tholatn=-1.
1042	tholats=0.
1043	tholone=0.
1044	tholonw=-1.
1045	alb_tho_spe=0.1 ! default value
1046	emis_tho_spe=1. ! default value
1047	call getin_p(" tholatn",tholatn)
1048	call getin_p(" tholats",tholats)
1049	call getin_p(" tholone",tholone)
1050	call getin_p(" tholonw",tholonw)
1051	if (is_master)write(,)trim(rname)//&
1052	" Values for special tholins albedo = ",tholatn,&
1053	tholats,tholone,tholonw,alb_tho_spe,emis_tho_spe
1054
1055	if (is_master)write(,)trim(rname)//&
1056	"Specific albedo"
1057	spelon1=-180. ! default value
1058	spelon2=180. ! default value
1059	spelat1=-90. ! default value
1060	spelat2=90. ! default value
1061	specalb=.false. ! default value
1062	if (is_master)write(,)trim(rname)//&
1063	"albedo/emis spe "
1064	albspe=0.1 ! default value
1065	emispe=1. ! default value
1066	call getin_p("spelon1",spelon1)
1067	call getin_p("spelon2",spelon2)
1068	call getin_p("spelat1",spelat1)
1069	call getin_p("spelat2",spelat2)
1070	call getin_p("specalb",specalb)
1071	call getin_p("albspe",albspe)
1072	call getin_p("emispe",emispe)
1073
1074	if (is_master)write(,)trim(rname)//&
1075	" specific = ",specalb
1076
1077	!******************** TI ***************************
1078
1079	if (is_master)write(,)trim(rname)//&
1080	"Change TI with time"
1081	changeti=.false. ! default value
1082	call getin_p("changeti",changeti)
1083	if (is_master)write(,)trim(rname)//&
1084	" changeti = ",changeti
1085	changetid=.false. ! default value for diurnal TI
1086	call getin_p("changetid",changetid)
1087	if (is_master)write(,)trim(rname)//&
1088	" changetid = ",changetid
1089
1090	if (is_master)write(,)trim(rname)//&
1091	"IT N2 "
1092	ITN2=800. ! default value
1093	call getin_p("ITN2",ITN2)
1094	if (is_master)write(,)trim(rname)//&
1095	" ITN2 = ",ITN2
1096	ITN2d=20. ! default value
1097	call getin_p("ITN2d",ITN2d)
1098	if (is_master)write(,)trim(rname)//&
1099	" ITN2d = ",ITN2d
1100
1101	if (is_master)write(,)trim(rname)//&
1102	"IT CH4"
1103	ITCH4=800. ! default value
1104	call getin_p("ITCH4",ITCH4)
1105	if (is_master)write(,)trim(rname)//&
1106	" ITCH4 = ",ITCH4
1107	ITCH4d=20. ! default value
1108	call getin_p("ITCH4d",ITCH4d)
1109	if (is_master)write(,)trim(rname)//&
1110	" ITCH4d = ",ITCH4d
1111
1112	if (is_master)write(,)trim(rname)//&
1113	"IT H2O"
1114	ITH2O=800. ! default value
1115	call getin_p("ITH2O",ITH2O)
1116	if (is_master)write(,)trim(rname)//&
1117	" ITH2O = ",ITH2O
1118	ITH2Od=20. ! default value
1119	call getin_p("ITH2Od",ITH2Od)
1120	if (is_master)write(,)trim(rname)//&
1121	" ITH2Od = ",ITH2Od
1122
1123	!************ COOLING *************
1124
1125	alpha_top=5e-11 ! default value
1126	call getin_p("alpha_top",alpha_top)
1127	if (is_master)write(,)trim(rname)//&
1128	" alpha_top = ",alpha_top
1129	pref=0.02 ! default value
1130	call getin_p("pref",pref)
1131	if (is_master)write(,)trim(rname)//&
1132	" pref = ",pref
1133	deltap=0.1 ! default value
1134	call getin_p("deltap",deltap)
1135	if (is_master)write(,)trim(rname)//&
1136	" deltap = ",deltap
1137
1138	!=================================
1139	! TWOLAY scheme and NH3 cloudare left for retrocompatibility only,
1140	! You should now use N-LAYER scheme (see below).
1141
1142	if (is_master) write(,)trim(rname)//&
1143	": Radiatively active two-layer aerosols?"
1144	aeroback2lay=.false. ! default value
1145	call getin_p("aeroback2lay",aeroback2lay)
1146	if (is_master) write(,)trim(rname)//": aeroback2lay = ",aeroback2lay
1147
1148	if (aeroback2lay.and.is_master) then
1149	print*,'Warning : The TWOLAY AEROSOL scheme is deprecated and buggy...'
1150	print*,'You should use the generic n-layer scheme (see aeronlay).'
1151	endif
1152
1153	if (is_master) write(,)trim(rname)//": Radiatively active ammonia cloud?"
1154	aeronh3=.false. ! default value
1155	call getin_p("aeronh3",aeronh3)
1156	if (is_master) write(,)trim(rname)//": aeronh3 = ",aeronh3
1157
1158	if (aeronh3.and.is_master) then
1159	print*,'Warning : You are using specific NH3 cloud scheme ...'
1160	print*,'You should use the generic n-layer scheme (see aeronlay).'
1161	endif
1162
1163	if (is_master) write(,)trim(rname)//&
1164	": TWOLAY AEROSOL: total optical depth "//&
1165	"in the tropospheric layer (visible)"
1166	obs_tau_col_tropo=8.D0
1167	call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
1168	if (is_master) write(,)trim(rname)//&
1169	": obs_tau_col_tropo = ",obs_tau_col_tropo
1170
1171	if (is_master) write(,)trim(rname)//&
1172	": TWOLAY AEROSOL: total optical depth "//&
1173	"in the stratospheric layer (visible)"
1174	obs_tau_col_strato=0.08D0
1175	call getin_p("obs_tau_col_strato",obs_tau_col_strato)
1176	if (is_master) write(,)trim(rname)//&
1177	": obs_tau_col_strato = ",obs_tau_col_strato
1178
1179	if (is_master) write(,)trim(rname)//&
1180	": TWOLAY AEROSOL: optprop_back2lay_vis?"
1181	optprop_back2lay_vis = 'optprop_saturn_vis_n20.dat'
1182	call getin_p("optprop_back2lay_vis",optprop_back2lay_vis)
1183	if (is_master) write(,)trim(rname)//&
1184	": optprop_back2lay_vis = ",trim(optprop_back2lay_vis)
1185
1186	if (is_master) write(,)trim(rname)//&
1187	": TWOLAY AEROSOL: optprop_back2lay_ir?"
1188	optprop_back2lay_ir = 'optprop_saturn_ir_n20.dat'
1189	call getin_p("optprop_back2lay_ir",optprop_back2lay_ir)
1190	if (is_master) write(,)trim(rname)//&
1191	": optprop_back2lay_ir = ",trim(optprop_back2lay_ir)
1192
1193	if (is_master) write(,)trim(rname)//&
1194	": TWOLAY AEROSOL: pres_bottom_tropo? in pa"
1195	pres_bottom_tropo=66000.0
1196	call getin_p("pres_bottom_tropo",pres_bottom_tropo)
1197	if (is_master) write(,)trim(rname)//&
1198	": pres_bottom_tropo = ",pres_bottom_tropo
1199
1200	if (is_master) write(,)trim(rname)//&
1201	": TWOLAY AEROSOL: pres_top_tropo? in pa"
1202	pres_top_tropo=18000.0
1203	call getin_p("pres_top_tropo",pres_top_tropo)
1204	if (is_master) write(,)trim(rname)//&
1205	": pres_top_tropo = ",pres_top_tropo
1206
1207	if (is_master) write(,)trim(rname)//&
1208	": TWOLAY AEROSOL: pres_bottom_strato? in pa"
1209	pres_bottom_strato=2000.0
1210	call getin_p("pres_bottom_strato",pres_bottom_strato)
1211	if (is_master) write(,)trim(rname)//&
1212	": pres_bottom_strato = ",pres_bottom_strato
1213
1214	! Sanity check
1215	if (pres_bottom_strato .gt. pres_top_tropo) then
1216	if(is_master) then
1217	print*,'Error : TWOLAY AEROSOL, Please ensure that in callphys.def'
1218	print*,'you have pres_top_tropo > pres_bottom_strato !'
1219	endif
1220	call abort_physic(rname," pres_top_tropo > pres_bottom_strato!",1)
1221	endif
1222
1223	if (is_master) write(,)trim(rname)//&
1224	": TWOLAY AEROSOL: pres_top_strato? in pa"
1225	pres_top_strato=100.0
1226	call getin_p("pres_top_strato",pres_top_strato)
1227	if (is_master) write(,)trim(rname)//&
1228	": pres_top_strato = ",pres_top_strato
1229
1230	if (is_master) write(,)trim(rname)//&
1231	": TWOLAY AEROSOL: particle size in the ", &
1232	"tropospheric layer, in meters"
1233	size_tropo=2.e-6
1234	call getin_p("size_tropo",size_tropo)
1235	if (is_master) write(,)trim(rname)//": size_tropo = ",size_tropo
1236
1237	if (is_master) write(,)trim(rname)//&
1238	": TWOLAY AEROSOL: particle size in the ", &
1239	"stratospheric layer, in meters"
1240	size_strato=1.e-7
1241	call getin_p("size_strato",size_strato)
1242	if (is_master) write(,)trim(rname)//": size_strato = ",size_strato
1243
1244	if (is_master) write(,)trim(rname)//&
1245	": NH3 (thin) cloud: total optical depth"
1246	tau_nh3_cloud=7.D0
1247	call getin_p("tau_nh3_cloud",tau_nh3_cloud)
1248	if (is_master) write(,)trim(rname)//": tau_nh3_cloud = ",tau_nh3_cloud
1249
1250	if (is_master) write(,)trim(rname)//": NH3 (thin) cloud pressure level"
1251	pres_nh3_cloud=7.D0
1252	call getin_p("pres_nh3_cloud",pres_nh3_cloud)
1253	if (is_master) write(,)trim(rname)//": pres_nh3_cloud = ",pres_nh3_cloud
1254
1255	if (is_master) write(,)trim(rname)//": NH3 (thin) cloud: particle sizes"
1256	size_nh3_cloud=3.e-6
1257	call getin_p("size_nh3_cloud",size_nh3_cloud)
1258	if (is_master) write(,)trim(rname)//": size_nh3_cloud = ",size_nh3_cloud
1259
1260	!=================================
1261	! Generic N-LAYER aerosol scheme
1262
1263	if (is_master) write(,)trim(rname)//&
1264	": Radiatively active generic n-layer aerosols?"
1265	aeronlay=.false. ! default value
1266	call getin_p("aeronlay",aeronlay)
1267	if (is_master) write(,)trim(rname)//": aeronlay = ",aeronlay
1268
1269	if (is_master) write(,)trim(rname)//&
1270	": Number of generic aerosols layers?"
1271	nlayaero=1 ! default value
1272	call getin_p("nlayaero",nlayaero)
1273	! Avoid to allocate arrays of size 0
1274	if (aeronlay .and. nlayaero.lt.1) then
1275	if (is_master) then
1276	print*, " You are trying to set no generic aerosols..."
1277	print*, " Set aeronlay=.false. instead ! I abort."
1278	endif
1279	call abort_physic(rname,"no generic aerosols...",1)
1280	endif
1281	if (.not. aeronlay) nlayaero=1
1282	if (is_master) write(,)trim(rname)//": nlayaero = ",nlayaero
1283
1284	! This is necessary, we just set the number of aerosol layers
1285	IF(.NOT.ALLOCATED(aeronlay_tauref)) ALLOCATE(aeronlay_tauref(nlayaero))
1286	IF(.NOT.ALLOCATED(aeronlay_lamref)) ALLOCATE(aeronlay_lamref(nlayaero))
1287	IF(.NOT.ALLOCATED(aeronlay_choice)) ALLOCATE(aeronlay_choice(nlayaero))
1288	IF(.NOT.ALLOCATED(aeronlay_pbot)) ALLOCATE(aeronlay_pbot(nlayaero))
1289	IF(.NOT.ALLOCATED(aeronlay_ptop)) ALLOCATE(aeronlay_ptop(nlayaero))
1290	IF(.NOT.ALLOCATED(aeronlay_sclhght)) ALLOCATE(aeronlay_sclhght(nlayaero))
1291	IF(.NOT.ALLOCATED(aeronlay_size)) ALLOCATE(aeronlay_size(nlayaero))
1292	IF(.NOT.ALLOCATED(aeronlay_nueff)) ALLOCATE(aeronlay_nueff(nlayaero))
1293	IF(.NOT.ALLOCATED(optprop_aeronlay_ir)) ALLOCATE(optprop_aeronlay_ir(nlayaero))
1294	IF(.NOT.ALLOCATED(optprop_aeronlay_vis)) ALLOCATE(optprop_aeronlay_vis(nlayaero))
1295
1296	if (is_master) write(,)trim(rname)//&
1297	": Generic n-layer aerosols: Optical depth at reference wavelenght"
1298	aeronlay_tauref=1.0E-1
1299	call getin_p("aeronlay_tauref",aeronlay_tauref)
1300	if (is_master) write(,)trim(rname)//&
1301	": aeronlay_tauref = ",aeronlay_tauref
1302
1303	if (is_master) write(,)trim(rname)//&
1304	": Generic n-layer aerosols: Reference wavelenght for optical depths (m)"
1305	aeronlay_lamref=0.6E-6
1306	call getin_p("aeronlay_lamref",aeronlay_lamref)
1307	if (is_master) write(,)trim(rname)//&
1308	": aeronlay_lamref = ",aeronlay_lamref
1309
1310	if (is_master) then
1311	write(,)trim(rname)//&
1312	": Generic n-layer aerosols: Vertical profile choice : "
1313	write(,)trim(rname)//&
1314	" (1) Tau btwn ptop and pbot follows atm. scale height"
1315	write(,)trim(rname)//&
1316	" or (2) Tau above pbot follows its own scale height"
1317	endif
1318	aeronlay_choice=1
1319	call getin_p("aeronlay_choice",aeronlay_choice)
1320	if (is_master) write(,)trim(rname)//&
1321	": aeronlay_choice = ",aeronlay_choice
1322
1323	if (is_master) write(,)trim(rname)//&
1324	": Generic n-layer aerosols: bottom pressures (Pa)"
1325	aeronlay_pbot=2000.0
1326	call getin_p("aeronlay_pbot",aeronlay_pbot)
1327	if (is_master) write(,)trim(rname)//": aeronlay_pbot = ",aeronlay_pbot
1328
1329	if (is_master) write(,)trim(rname)//&
1330	": Generic n-layer aerosols: (if choice=1) Top pressures (Pa) "
1331	aeronlay_ptop=300000.0
1332	call getin_p("aeronlay_ptop",aeronlay_ptop)
1333	if (is_master) write(,)trim(rname)//": aeronlay_ptop = ",aeronlay_ptop
1334
1335	if (is_master) write(,)trim(rname)//&
1336	": Generic n-layer aerosols: (if choice=2) Scale height / atm. scale height"
1337	aeronlay_sclhght=0.2
1338	call getin_p("aeronlay_sclhght",aeronlay_sclhght)
1339	if (is_master) write(,)trim(rname)//&
1340	": aeronlay_sclhght = ",aeronlay_sclhght
1341
1342	if (is_master) write(,)trim(rname)//&
1343	": Generic n-layer aerosols: particles effective radii (m)"
1344	aeronlay_size=1.e-6
1345	call getin_p("aeronlay_size",aeronlay_size)
1346	if (is_master) write(,)trim(rname)//": aeronlay_size = ",aeronlay_size
1347
1348	if (is_master) write(,)trim(rname)//&
1349	": Generic n-layer aerosols: particles radii effective variance"
1350	aeronlay_nueff=0.1
1351	call getin_p("aeronlay_nueff",aeronlay_nueff)
1352	if (is_master) write(,)trim(rname)//": aeronlay_nueff = ",aeronlay_nueff
1353
1354	if (is_master) write(,)trim(rname)//&
1355	": Generic n-layer aerosols: VIS optical properties file"
1356	optprop_aeronlay_vis = 'optprop_saturn_vis_n20.dat'
1357	call getin_p("optprop_aeronlay_vis",optprop_aeronlay_vis)
1358	if (is_master) write(,)trim(rname)//&
1359	": optprop_aeronlay_vis = ",optprop_aeronlay_vis
1360
1361	if (is_master) write(,)trim(rname)//&
1362	": Generic n-layer aerosols: IR optical properties file"
1363	optprop_aeronlay_ir = 'optprop_saturn_ir_n20.dat'
1364	call getin_p("optprop_aeronlay_ir",optprop_aeronlay_ir)
1365	if (is_master) write(,)trim(rname)//&
1366	": optprop_aeronlay_ir = ",optprop_aeronlay_ir
1367
1368
1369	!=================================
1370
1371	if (is_master) write(,)trim(rname)//&
1372	": Is the variable gas species radiatively active?"
1373	Tstrat=167.0
1374	varactive=.false.
1375	call getin_p("varactive",varactive)
1376	if (is_master) write(,)trim(rname)//": varactive = ",varactive
1377
1378	if (is_master) write(,)trim(rname)//&
1379	": Is the variable gas species distribution set?"
1380	varfixed=.false.
1381	call getin_p("varfixed",varfixed)
1382	if (is_master) write(,)trim(rname)//": varfixed = ",varfixed
1383
1384	if (is_master) write(,)trim(rname)//&
1385	": What is the saturation % of the variable species?"
1386	satval=0.8
1387	call getin_p("satval",satval)
1388	if (is_master) write(,)trim(rname)//": satval = ",satval
1389
1390
1391	! Test of incompatibility:
1392	! if varactive, then varfixed should be false
1393	if (varactive.and.varfixed) then
1394	call abort_physic(rname,'if varactive, varfixed must be OFF!',1)
1395	endif
1396
1397	if (is_master) write(,)trim(rname)//": Gravitationnal sedimentation ?"
1398	sedimentation=.false. ! default value
1399	call getin_p("sedimentation",sedimentation)
1400	if (is_master) write(,)trim(rname)//": sedimentation = ",sedimentation
1401
1402	if (is_master) write(,)trim(rname)//": Generic Condensation of tracers ?"
1403	generic_condensation=.false. !default value
1404	call getin_p("generic_condensation",generic_condensation)
1405	if (is_master) write(,)trim(rname)//": generic_condensation = ",generic_condensation
1406
1407	if (is_master) write(,)trim(rname)//": Spectral Dependant albedo ?"
1408	albedo_spectral_mode=.false. ! default value
1409	call getin_p("albedo_spectral_mode",albedo_spectral_mode)
1410	if (is_master) write(,)trim(rname)//&
1411	": albedo_spectral_mode = ",albedo_spectral_mode
1412
1413	if (is_master) then
1414	write(,)trim(rname)//": Snow albedo ?"
1415	write(,)trim(rname)//": If albedo_spectral_mode=.true., then this "
1416	write(,)trim(rname)//": corresponds to the 0.5 microns snow albedo."
1417	endif
1418	albedosnow=0.5 ! default value
1419	call getin_p("albedosnow",albedosnow)
1420	if (is_master) write(,)trim(rname)//": albedosnow = ",albedosnow
1421
1422	if (is_master) write(,)trim(rname)//": N2 ice albedo ?"
1423	albedon2ice=0.5 ! default value
1424	call getin_p("albedon2ice",albedon2ice)
1425	if (is_master) write(,)trim(rname)//": albedon2ice = ",albedon2ice
1426
1427	if (is_master) write(,)trim(rname)//": Maximum ice thickness ?"
1428	maxicethick=2.0 ! default value
1429	call getin_p("maxicethick",maxicethick)
1430	if (is_master) write(,)trim(rname)//": maxicethick = ",maxicethick
1431
1432	if (is_master) write(,)trim(rname)//": Minimum eddy mix coeff in 1D ?"
1433	kmixmin=1.0e-2 ! default value
1434	call getin_p("kmixmin",kmixmin)
1435	if (is_master) write(,)trim(rname)//": kmixmin = ",kmixmin
1436
1437	if (is_master) write(,)'Predefined Cp from dynamics is ',cpp,'J kg^-1 K^-1'
1438	if (is_master) write(,)'Predefined Mg from dynamics is ',mugaz,'amu'
1439
1440	force_cpp=.false. ! default value
1441	call getin_p("force_cpp",force_cpp)
1442	if (force_cpp) then
1443	if (is_master) write(,)trim(rname)//": force_cpp = ",force_cpp
1444	if (is_master) write(,)trim(rname)//": force_cpp is deprecated.",&
1445	"Set cpp_mugaz_mode=1 in callfis to emulate force_cpp=.true."
1446	call abort_physic(rname,"Anyway, you need to set force_cpp=.false. to continue.",1)
1447	endif
1448
1449	if (is_master) write(,)trim(rname)//&
1450	": where do you want your cpp/mugaz value to come from?",&
1451	"=> 0: dynamics (3d), 1: forced in callfis (1d), 2: computed from gases.def (1d)?"
1452	cpp_mugaz_mode = 0 ! default value
1453	call getin_p("cpp_mugaz_mode",cpp_mugaz_mode)
1454	if (is_master) write(,)trim(rname)//": cpp_mugaz_mode = ",cpp_mugaz_mode
1455
1456
1457	! Test of incompatibility:
1458
1459	if ((.not.tracer).and.(haze)) then
1460	call abort_physic(rname, 'if haze are on, tracers must be on!', 1)
1461	endif
1462	if (haze_proffix.and.sedimentation) then
1463	call abort_physic(rname, 'if haze profile is set, sedimentation must be deactivated', 1)
1464	endif
1465	if (callmolvis.and..not.callconduct) then
1466	call abort_physic(rname, 'if callmolvis is set, callconduct must be true', 1)
1467	endif
1468	if (glaflow.and..not.fast) then
1469	call abort_physic(rname, 'if glaflow is set, fast must be true', 1)
1470	endif
1471	if (paleo.and..not.fast) then
1472	call abort_physic(rname, 'if paleo is set, fast must be true', 1)
1473	endif
1474	if ((haze_proffix.or.haze_radproffix).and..not.aerohaze) then
1475	call abort_physic(rname, 'for now, haze/rad proffix only works w aerohaze=T', 1)
1476	endif
1477	if (carbox.and.condcosurf.and.no_n2frost) then
1478	call abort_physic(rname, "CO surface condensation and no_n2frost are both active which may not be relevant", 1)
1479	end if
1480
1481	if ((cpp_mugaz_mode >= 1).and.(is_master).and.(ngrid>1)) then
1482	write(,)' !!!!! Be aware that having different values for cpp and mugaz in the dynamics and physics'
1483	write(,)' in 3D can result in very pathological behavior. You have been warned !!!!!'
1484	else if ((cpp_mugaz_mode == 0).and.(is_master).and.(ngrid==1)) then
1485	! for this specific 1D error we will remove run.def before aborting JL22
1486	call system("rm -rf run.def")
1487	call abort_physic(rname,"cpp_mugaz_mode must be >= 1 in 1d",1)
1488	endif
1489
1490	if (cpp_mugaz_mode == 1) then
1491	mugaz = -99999.
1492	if (is_master) write(,)trim(rname)//&
1493	": MEAN MOLECULAR MASS in g mol-1 ?"
1494	call getin_p("mugaz",mugaz)
1495	IF (mugaz.eq.-99999.) THEN
1496	call abort_physic(rname,"mugaz must be set if cpp_mugaz_mode = 1",1)
1497	ENDIF
1498	cpp = -99999.
1499	if (is_master) write(,)trim(rname)//&
1500	": SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
1501	call getin_p("cpp",cpp)
1502	IF (cpp.eq.-99999.) THEN
1503	PRINT *, "cpp must be set if cpp_mugaz_mode = 1"
1504	STOP
1505	ENDIF
1506	if (is_master) write(,)'New Cp from callfis is ',cpp,'J kg^-1 K^-1'
1507	if (is_master) write(,)'New Mg from callfis is ',mugaz,'amu'
1508
1509	endif ! of if (cpp_mugaz_mode == 1)
1510	call su_gases
1511	call calc_cpp_mugaz
1512
1513	if (is_master) then
1514	PRINT*,'--------------------------------------------'
1515	PRINT*
1516	PRINT*
1517	endif
1518	ELSE
1519	call abort_physic(rname,'Cannot read file callphys.def. Is it here ?',1)
1520	ENDIF ! of IF(iscallphys)
1521
1522	if (is_master) then
1523	PRINT*
1524	PRINT*,'inifis: daysec',daysec
1525	PRINT*
1526	PRINT*,'inifis: The radiative transfer is computed:'
1527	PRINT*,' each ',iradia,' physical time-step'
1528	PRINT,' or each ',iradiadtphys,' seconds'
1529	PRINT*
1530	endif
1531
1532	!-----------------------------------------------------------------------
1533	! Some more initialization:
1534	! ------------------------
1535
1536	! Initializations for comgeomfi_h
1537	#ifndef MESOSCALE
1538	totarea=SSUM(ngrid,parea,1)
1539	call planetwide_sumval(parea,totarea_planet)
1540
1541	!! those are defined in comdiurn_h.F90
1542	IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
1543	IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
1544	IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
1545	IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
1546
1547	DO ig=1,ngrid
1548	sinlat(ig)=sin(plat(ig))
1549	coslat(ig)=cos(plat(ig))
1550	sinlon(ig)=sin(plon(ig))
1551	coslon(ig)=cos(plon(ig))
1552	ENDDO
1553	#endif
1554	! initialize variables in radinc_h
1555	call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
1556
1557	! initialize variables and allocate arrays in radcommon_h
1558	call ini_radcommon_h(naerkind)
1559
1560	! allocate "comsoil_h" arrays
1561	call ini_comsoil_h(ngrid)
1562
1563	END SUBROUTINE inifis
1564
1565	END MODULE inifis_mod

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

Download in other formats:

Original Format