Context Navigation

inifis_mod.F90 @ 3453

Last change on this file since 3453 was 3439, checked in by afalco, 9 months ago
Pluto PCM: added some incompatibility checks. AF
File size: 56.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	! Slab Ocean !AF24: removed
409
410	! Chemistry in the thermosphere
411	if (is_master) write(,) trim(rname)//": Use deposition ?"
412	depos=.false. ! default value
413	call getin_p("depos",depos)
414	if (is_master) write(,) trim(rname)//": depos = ",depos
415
416	if (is_master) write(,)trim(rname)//": Production of haze ?"
417	haze=.false. ! default value
418	call getin_p("haze",haze)
419	if (is_master) write(,)trim(rname)//": haze = ",haze
420
421	! Global1D mean and solar zenith angle
422
423	if (ngrid.eq.1) then
424	PRINT*, 'Simulate global averaged conditions ?'
425	global1d = .false. ! default value
426	call getin_p("global1d",global1d)
427	write(,) "global1d = ",global1d
428
429	! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
430	if (global1d.and.diurnal) then
431	call abort_physic(rname,'if global1d is true, diurnal must be set to false',1)
432	endif
433
434	if (global1d) then
435	PRINT *,'Solar Zenith angle (deg.) ?'
436	PRINT *,'(assumed for averaged solar flux S/4)'
437	szangle=60.0 ! default value
438	call getin_p("szangle",szangle)
439	write(,) "szangle = ",szangle
440	endif
441	endif ! of if (ngrid.eq.1)
442
443	! Test of incompatibility:
444	! if kastprof used, we must be in 1D
445	if (kastprof.and.(ngrid.gt.1)) then
446	call abort_physic(rname,'kastprof can only be used in 1D!',1)
447	endif
448
449	if (is_master) write(,)trim(rname)//&
450	": Stratospheric temperature for kastprof mode?"
451	Tstrat=167.0
452	call getin_p("Tstrat",Tstrat)
453	if (is_master) write(,)trim(rname)//": Tstrat = ",Tstrat
454
455	if (is_master) write(,)trim(rname)//": Remove lower boundary?"
456	nosurf=.false.
457	call getin_p("nosurf",nosurf)
458	if (is_master) write(,)trim(rname)//": nosurf = ",nosurf
459
460	! Tests of incompatibility:
461	if (nosurf.and.callsoil) then
462	if (is_master) then
463	write(,)trim(rname)//'nosurf not compatible with soil scheme!'
464	write(,)trim(rname)//'... got to make a choice!'
465	endif
466	call abort_physic(rname,"nosurf not compatible with soil scheme!",1)
467	endif
468
469	if (is_master) write(,)trim(rname)//": Add an internal heat flux?", &
470	"... matters only if callsoil=F"
471	intheat=0.
472	call getin_p("intheat",intheat)
473	if (is_master) write(,)trim(rname)//": intheat = ",intheat
474
475
476	if (is_master) write(,)trim(rname)//": Test physics timescale in 1D?"
477	testradtimes=.false.
478	call getin_p("testradtimes",testradtimes)
479	if (is_master) write(,)trim(rname)//": testradtimes = ",testradtimes
480
481	! Test of incompatibility:
482	! if testradtimes used, we must be in 1D
483	if (testradtimes.and.(ngrid.gt.1)) then
484	call abort_physic(rname,'testradtimes can only be used in 1D!',1)
485	endif
486
487	if (is_master) write(,)trim(rname)//": Default planetary temperature?"
488	tplanet=215.0
489	call getin_p("tplanet",tplanet)
490	if (is_master) write(,)trim(rname)//": tplanet = ",tplanet
491
492	if (is_master) write(,)trim(rname)//": Which star?"
493	startype=1 ! default value = Sol
494	call getin_p("startype",startype)
495	if (is_master) write(,)trim(rname)//": startype = ",startype
496
497	if (is_master) write(,)trim(rname)//": Value of stellar flux at 1 AU?"
498	Fat1AU=1356.0 ! default value = Sol today
499	call getin_p("Fat1AU",Fat1AU)
500	if (is_master) write(,)trim(rname)//": Fat1AU = ",Fat1AU
501
502
503	! TRACERS:
504
505	if (is_master) write(,)trim(rname)//": fixed radii for Cloud particles?"
506	radfixed=.false. ! default value
507	call getin_p("radfixed",radfixed)
508	if (is_master) write(,)trim(rname)//": radfixed = ",radfixed
509
510	if(kastprof)then
511	radfixed=.true.
512	endif
513
514	if (is_master) write(,)trim(rname)//&
515	"Number of radiatively active aerosols:"
516	naerkind=0 ! default value
517	call getin_p("naerkind",naerkind)
518	if (is_master) write(,)trim(rname)//": naerkind = ",naerkind
519
520
521	!***************************************************************
522	!! TRACERS options
523
524	if (is_master)write(,)trim(rname)//&
525	"call N2 condensation ?"
526	n2cond=.true. ! default value
527	call getin_p("n2cond",n2cond)
528	if (is_master)write(,)trim(rname)//&
529	" n2cond = ",n2cond
530
531	if (is_master)write(,)trim(rname)//&
532	"call N2 cloud condensation ?"
533	condensn2=.false. ! default value
534	call getin_p("condensn2",condensn2)
535	if (is_master)write(,)trim(rname)//&
536	"condensn2 = ",condensn2
537
538	if (is_master)write(,)trim(rname)//&
539	"call no N2 frost formation?"
540	no_n2frost=.false. ! default value
541	call getin_p("no_n2frost",no_n2frost)
542	if (is_master)write(,)trim(rname)//&
543	"no_n2frost = ",no_n2frost
544
545	if (is_master)write(,)trim(rname)//&
546	"N2 condensation subtimestep?"
547	nbsub=20 ! default value
548	call getin_p("nbsub",nbsub)
549	if (is_master)write(,)trim(rname)//&
550	" nbsub = ",nbsub
551
552	if (is_master)write(,)trim(rname)//&
553	"Gravitationnal sedimentation ?"
554	sedimentation=.true. ! default value
555	call getin_p("sedimentation",sedimentation)
556	if (is_master)write(,)trim(rname)//&
557	" sedimentation = ",sedimentation
558
559	if (is_master)write(,)trim(rname)//&
560	"Compute glacial flow ?"
561	glaflow=.false. ! default value
562	call getin_p("glaflow",glaflow)
563	if (is_master)write(,)trim(rname)//&
564	" glaflow = ",glaflow
565
566	if (is_master)write(,)trim(rname)//&
567	"Compute methane cycle ?"
568	methane=.false. ! default value
569	call getin_p("methane",methane)
570	if (is_master)write(,)trim(rname)//&
571	" methane = ",methane
572	condmetsurf=.true. ! default value
573	call getin_p("condmetsurf",condmetsurf)
574	if (is_master)write(,)trim(rname)//&
575	" condmetsurf = ",condmetsurf
576
577	if (is_master)write(,)trim(rname)//&
578	"Compute methane clouds ?"
579	metcloud=.false. ! default value
580	call getin_p("metcloud",metcloud)
581	if (is_master)write(,)trim(rname)//&
582	" metcloud = ",metcloud
583
584	if (is_master)write(,)trim(rname)//&
585	"Compute CO cycle ?"
586	carbox=.false. ! default value
587	call getin_p("carbox",carbox)
588	if (is_master)write(,)trim(rname)//&
589	" carbox = ",carbox
590	condcosurf=.true. ! default value
591	call getin_p("condcosurf",condcosurf)
592	if (is_master)write(,)trim(rname)//&
593	" condcosurf = ",condcosurf
594
595	if (is_master)write(,)trim(rname)//&
596	"Compute CO clouds ?"
597	monoxcloud=.false. ! default value
598	call getin_p("monoxcloud",monoxcloud)
599	if (is_master)write(,)trim(rname)//&
600	" monoxcloud = ",monoxcloud
601
602	if (is_master)write(,)trim(rname)//&
603	"atmospheric redistribution (s):"
604	tau_n2=1. ! default value
605	call getin_p("tau_n2",tau_n2)
606	if (is_master)write(,)trim(rname)//&
607	" tau_n2 = ",tau_n2
608	tau_ch4=1.E7 ! default value
609	call getin_p("tau_ch4",tau_ch4)
610	if (is_master)write(,)trim(rname)//&
611	" tau_ch4 = ",tau_ch4
612	tau_co=1. ! default value
613	call getin_p("tau_co",tau_co)
614	if (is_master)write(,)trim(rname)//&
615	" tau_co = ",tau_co
616	tau_prechaze=1. ! default value
617	call getin_p("tau_prechaze",tau_prechaze)
618	if (is_master)write(,)trim(rname)//&
619	" tau_prechaze = ",tau_prechaze
620
621	if (is_master)write(,)trim(rname)//&
622	"Day fraction for limited cold trap in SP?"
623	dayfrac=0. ! default value
624	call getin_p("dayfrac",dayfrac)
625	if (is_master)write(,)trim(rname)//&
626	" dayfrac = ",dayfrac
627	thresh_non2=0. ! default value
628	call getin_p("thresh_non2",thresh_non2)
629	if (is_master)write(,)trim(rname)//&
630	" thresh_non2 = ",thresh_non2
631
632	! use Pluto.old routines
633	if (is_master) write(,) trim(rname)//": use vdifc from old Pluto ?"
634	oldplutovdifc=.false. ! default value
635	call getin_p("oldplutovdifc",oldplutovdifc)
636	if (is_master) write(,) trim(rname)//": oldplutovdifc = ",oldplutovdifc
637
638	if (is_master) write(,) trim(rname)//&
639	": call pluto.old correlated-k radiative transfer ?"
640	oldplutocorrk=.false. ! default value
641	call getin_p("oldplutocorrk",oldplutocorrk)
642	if (is_master) write(,) trim(rname)//": oldplutocorrk = ",oldplutocorrk
643
644	if (is_master) write(,) trim(rname)//&
645	": call pluto.old sedimentation ?"
646	oldplutosedim=.false. ! default value
647	call getin_p("oldplutosedim",oldplutosedim)
648	if (is_master) write(,) trim(rname)//": oldplutosedim = ",oldplutosedim
649
650	!***************************************************************
651	!! Haze options
652
653	if (is_master)write(,)trim(rname)//&
654	"Production of haze ?"
655	haze=.false. ! default value
656	call getin_p("haze",haze)
657	if (is_master)write(,)trim(rname)//&
658	" haze = ",haze
659	hazeconservch4=.false. ! conservch4, by default value ch4 is photolyzed
660	call getin_p("hazeconservch4",hazeconservch4)
661	if (is_master)write(,)trim(rname)//&
662	"hazconservch4 = ",hazeconservch4
663	if (is_master)write(,)trim(rname)//&
664	"Production of haze (fast version) ?"
665	fasthaze=.false. ! default value
666	call getin_p("fasthaze",fasthaze)
667	if (is_master)write(,)trim(rname)//&
668	"fasthaze = ",fasthaze
669
670	if (is_master)write(,)trim(rname)//&
671	"Add source of haze ?"
672	source_haze=.false. ! default value
673	call getin_p("source_haze",source_haze)
674	if (is_master)write(,)trim(rname)//&
675	" source_haze = ",source_haze
676	mode_hs=0 ! mode haze source default value
677	call getin_p("mode_hs",mode_hs)
678	if (is_master)write(,)trim(rname)//&
679	" mode_hs",mode_hs
680	kfix=1 ! default value
681	call getin_p("kfix",kfix)
682	if (is_master)write(,)trim(rname)//&
683	"levels kfix",kfix
684	fracsource=0.1 ! default value
685	call getin_p("fracsource",fracsource)
686	if (is_master)write(,)trim(rname)//&
687	" fracsource",fracsource
688	latsource=30. ! default value
689	call getin_p("latsource",latsource)
690	if (is_master)write(,)trim(rname)//&
691	" latsource",latsource
692	lonsource=180. ! default value
693	call getin_p("lonsource",lonsource)
694	if (is_master)write(,)trim(rname)//&
695	" lonsource",lonsource
696
697	if (is_master)write(,)trim(rname)//&
698	"Radiatively active haze ?"
699	aerohaze=.false. ! default value
700	call getin_p("aerohaze",aerohaze)
701	if (is_master)write(,)trim(rname)//&
702	"aerohaze = ",aerohaze
703
704	if (is_master)write(,)trim(rname)//&
705	"Haze monomer radius ?"
706	rad_haze=10.e-9 ! default value
707	call getin_p("rad_haze",rad_haze)
708	if (is_master)write(,)trim(rname)//&
709	"rad_haze = ",rad_haze
710
711	if (is_master)write(,)trim(rname)//&
712	"fractal particle ?"
713	fractal=.false. ! default value
714	call getin_p("fractal",fractal)
715	if (is_master)write(,)trim(rname)//&
716	"fractal = ",fractal
717	nb_monomer=10 ! default value
718	call getin_p("nb_monomer",nb_monomer)
719	if (is_master)write(,)trim(rname)//&
720	"nb_monomer = ",nb_monomer
721
722	if (is_master)write(,)trim(rname)//&
723	"fixed gaseous CH4 mixing ratio for rad transfer ?"
724	ch4fix=.false. ! default value
725	call getin_p("ch4fix",ch4fix)
726	if (is_master)write(,)trim(rname)//&
727	" ch4fix = ",ch4fix
728	if (is_master)write(,)trim(rname)//&
729	"fixed gaseous CH4 mixing ratio for rad transfer ?"
730	vmrch4fix=0.5 ! default value
731	call getin_p("vmrch4fix",vmrch4fix)
732	if (is_master)write(,)trim(rname)//&
733	" vmrch4fix = ",vmrch4fix
734	vmrch4_proffix=.false. ! default value
735	call getin_p("vmrch4_proffix",vmrch4_proffix)
736	if (is_master)write(,)trim(rname)//&
737	" vmrch4_proffix = ",vmrch4_proffix
738
739	if (is_master)write(,)trim(rname)//&
740	"call specific cooling for radiative transfer ?"
741	cooling=.false. ! default value
742	call getin_p("cooling",cooling)
743	if (is_master)write(,)trim(rname)//&
744	" cooling = ",cooling
745
746	if (is_master)write(,)trim(rname)//&
747	"NLTE correction for methane heating rates?"
748	nlte=.false. ! default value
749	call getin_p("nlte",nlte)
750	if (is_master)write(,)trim(rname)//&
751	" nlte = ",nlte
752	strobel=.false. ! default value
753	call getin_p("strobel",strobel)
754	if (is_master)write(,)trim(rname)//&
755	" strobel = ",strobel
756
757	if (is_master)write(,)trim(rname)//&
758	"fixed radius profile from txt file ?"
759	haze_radproffix=.false. ! default value
760	call getin_p("haze_radproffix",haze_radproffix)
761	if (is_master)write(,)trim(rname)//&
762	"haze_radproffix = ",haze_radproffix
763	if (is_master)write(,)trim(rname)//&
764	"fixed MMR profile from txt file ?"
765	haze_proffix=.false. ! default value
766	call getin_p("haze_proffix",haze_proffix)
767	if (is_master)write(,)trim(rname)//&
768	"haze_proffix = ",haze_proffix
769
770	if (is_master)write(,)trim(rname)//&
771	"Number mix ratio of haze particles for co clouds:"
772	Nmix_co=100000. ! default value
773	call getin_p("Nmix_co",Nmix_co)
774	if (is_master)write(,)trim(rname)//&
775	" Nmix_co = ",Nmix_co
776
777	if (is_master)write(,)trim(rname)//&
778	"Number mix ratio of haze particles for ch4 clouds:"
779	Nmix_ch4=100000. ! default value
780	call getin_p("Nmix_ch4",Nmix_ch4)
781	if (is_master)write(,)trim(rname)//&
782	" Nmix_ch4 = ",Nmix_ch4
783
784	!***************************************************************
785	!! Surface properties
786
787	!********* N2 *******************************
788
789	if (is_master)write(,)trim(rname)//&
790	"Mode for changing N2 albedo"
791	mode_n2=0 ! default value
792	call getin_p("mode_n2",mode_n2)
793	if (is_master)write(,)trim(rname)//&
794	" mode_n2 = ",mode_n2
795	thres_n2ice=1. ! default value
796	call getin_p("thres_n2ice",thres_n2ice)
797	if (is_master)write(,)trim(rname)//&
798	" thres_n2ice = ",thres_n2ice
799
800	if (is_master)write(,)trim(rname)//&
801	"Diff of N2 albedo with thickness"
802	deltab=0. ! default value
803	call getin_p("deltab",deltab)
804	if (is_master)write(,)trim(rname)//&
805	" deltab = ",deltab
806
807	if (is_master)write(,)trim(rname)//&
808	"albedo N2 beta "
809	alb_n2b=0.7 ! default value
810	call getin_p("alb_n2b",alb_n2b)
811	if (is_master)write(,)trim(rname)//&
812	" alb_n2b = ",alb_n2b
813
814	if (is_master)write(,)trim(rname)//&
815	"albedo N2 alpha "
816	alb_n2a=0.7 ! default value
817	call getin_p("alb_n2a",alb_n2a)
818	if (is_master)write(,)trim(rname)//&
819	" alb_n2a = ",alb_n2a
820
821	if (is_master)write(,)trim(rname)//&
822	"emis N2 beta "
823	emis_n2b=0.7 ! default value
824	call getin_p("emis_n2b",emis_n2b)
825	if (is_master)write(,)trim(rname)//&
826	" emis_n2b = ",emis_n2b
827
828	if (is_master)write(,)trim(rname)//&
829	"emis N2 alpha "
830	emis_n2a=0.7 ! default value
831	call getin_p("emis_n2a",emis_n2a)
832	if (is_master)write(,)trim(rname)//&
833	" emis_n2a = ",emis_n2a
834
835	!********* CH4 *******************************
836
837	if (is_master)write(,)trim(rname)//&
838	"Mode for changing CH4 albedo"
839	mode_ch4=0 ! default value
840	call getin_p("mode_ch4",mode_ch4)
841	if (is_master)write(,)trim(rname)//&
842	" mode_ch4 = ",mode_ch4
843	feedback_met=0 ! default value
844	call getin_p("feedback_met",feedback_met)
845	if (is_master)write(,)trim(rname)//&
846	" feedback_met = ",feedback_met
847	thres_ch4ice=1. ! default value, mm
848	call getin_p("thres_ch4ice",thres_ch4ice)
849	if (is_master)write(,)trim(rname)//&
850	" thres_ch4ice = ",thres_ch4ice
851	fdch4_latn=-1.
852	fdch4_lats=0.
853	fdch4_lone=0.
854	fdch4_lonw=-1.
855	fdch4_depalb=0.5
856	fdch4_finalb=0.95
857	fdch4_maxalb=0.99
858	fdch4_ampl=200.
859	fdch4_maxice=100.
860	call getin_p("fdch4_latn",fdch4_latn)
861	call getin_p("fdch4_lats",fdch4_lats)
862	call getin_p("fdch4_lone",fdch4_lone)
863	call getin_p("fdch4_lonw",fdch4_lonw)
864	call getin_p("fdch4_depalb",fdch4_depalb)
865	call getin_p("fdch4_finalb",fdch4_finalb)
866	call getin_p("fdch4_maxalb",fdch4_maxalb)
867	call getin_p("fdch4_maxice",fdch4_maxice)
868	call getin_p("fdch4_ampl",fdch4_ampl)
869	if (is_master)write(,)trim(rname)//&
870	" Values for albedo feedback = ",fdch4_latn,&
871	fdch4_lats,fdch4_lone,fdch4_lonw,fdch4_depalb,&
872	fdch4_finalb,fdch4_maxalb,fdch4_maxice,fdch4_ampl
873
874	if (is_master)write(,)trim(rname)//&
875	"Latitude for diff albedo methane"
876	metlateq=25. ! default value
877	call getin_p("metlateq",metlateq)
878	if (is_master)write(,)trim(rname)//&
879	" metlateq = ",metlateq
880
881	if (is_master)write(,)trim(rname)//&
882	"Ls1 and Ls2 for change of ch4 albedo"
883	metls1=-1. ! default value
884	metls2=-2. ! default value
885	call getin_p("metls1",metls1)
886	call getin_p("metls2",metls2)
887
888	if (is_master)write(,)trim(rname)//&
889	"albedo CH4 "
890	alb_ch4=0.5 ! default value
891	call getin_p("alb_ch4",alb_ch4)
892	if (is_master)write(,)trim(rname)//&
893	" alb_ch4 = ",alb_ch4
894
895	if (is_master)write(,)trim(rname)//&
896	"albedo equatorial CH4 "
897	alb_ch4_eq=alb_ch4 ! default value
898	call getin_p("alb_ch4_eq",alb_ch4_eq)
899	if (is_master)write(,)trim(rname)//&
900	" alb_ch4_eq = ",alb_ch4_eq
901
902	if (is_master)write(,)trim(rname)//&
903	"albedo south hemis CH4 "
904	alb_ch4_s=alb_ch4 ! default value
905	call getin_p("alb_ch4_s",alb_ch4_s)
906	if (is_master)write(,)trim(rname)//&
907	" alb_ch4_s = ",alb_ch4_s
908
909	if (is_master)write(,)trim(rname)//&
910	"emis CH4 "
911	emis_ch4=0.5 ! default value
912	call getin_p("emis_ch4",emis_ch4)
913	if (is_master)write(,)trim(rname)//&
914	" emis_ch4 = ",emis_ch4
915
916	if (is_master)write(,)trim(rname)//&
917	"CH4 lag for n2 sublimation limitation"
918	ch4lag=.false. ! default value
919	latlag=-90. ! default value
920	vmrlag=1. ! default value
921	call getin_p("ch4lag",ch4lag)
922	call getin_p("latlag",latlag)
923	call getin_p("vmrlag",vmrlag)
924	if (is_master)write(,)trim(rname)//&
925	" ch4lag = ",ch4lag
926	if (is_master)write(,)trim(rname)//&
927	" latlag = ",latlag
928	if (is_master)write(,)trim(rname)//&
929	" vmrlag = ",vmrlag
930
931	if (is_master)write(,)trim(rname)//&
932	"Max temperature for surface ?"
933	tsurfmax=.false. ! default value
934	albmin_ch4=0.3 ! default value
935	call getin_p("tsurfmax",tsurfmax)
936	call getin_p("albmin_ch4",albmin_ch4)
937	if (is_master)write(,)trim(rname)//&
938	" tsurfmax = ",tsurfmax
939	if (is_master)write(,)trim(rname)//&
940	" albmin_ch4 = ",albmin_ch4
941
942
943	if (is_master)write(,)trim(rname)//&
944	"fixed gaseous CH4 mixing ratio for rad transfer ?"
945	ch4fix=.false. ! default value
946	call getin_p("ch4fix",ch4fix)
947	if (is_master)write(,)trim(rname)//&
948	" ch4fix = ",ch4fix
949	if (is_master)write(,)trim(rname)//&
950	"fixed gaseous CH4 mixing ratio for rad transfer ?"
951	vmrch4fix=0.5 ! default value
952	call getin_p("vmrch4fix",vmrch4fix)
953	if (is_master)write(,)trim(rname)//&
954	" vmrch4fix = ",vmrch4fix
955	vmrch4_proffix=.false. ! default value
956	call getin_p("vmrch4_proffix",vmrch4_proffix)
957	if (is_master)write(,)trim(rname)//&
958	" vmrch4_proffix = ",vmrch4_proffix
959
960
961	!********* CO *******************************
962
963	if (is_master)write(,)trim(rname)//&
964	"albedo CO "
965	alb_co=0.4 ! default value
966	call getin_p("alb_co",alb_co)
967	if (is_master)write(,)trim(rname)//&
968	" alb_co = ",alb_co
969	thres_coice=1. ! default value, mm
970	call getin_p("thres_coice",thres_coice)
971	if (is_master)write(,)trim(rname)//&
972	" thres_coice = ",thres_coice
973
974	if (is_master)write(,)trim(rname)//&
975	"emis CO "
976	emis_co=0.4 ! default value
977	call getin_p("emis_co",emis_co)
978	if (is_master)write(,)trim(rname)//&
979	" emis_co = ",emis_co
980
981	!********* THOLINS *******************************
982	if (is_master)write(,)trim(rname)//&
983	"Mode for changing tholins albedo/emis"
984	mode_tholins=0 ! default value
985	call getin_p("mode_tholins",mode_tholins)
986	if (is_master)write(,)trim(rname)//&
987	" mode_tholins = ",mode_tholins
988
989	if (is_master)write(,)trim(rname)//&
990	"albedo tho "
991	alb_tho=0.1 ! default value
992	call getin_p("alb_tho",alb_tho)
993	if (is_master)write(,)trim(rname)//&
994	" alb_tho = ",alb_tho
995
996	if (is_master)write(,)trim(rname)//&
997	"albedo tho eq"
998	alb_tho_eq=0.1 ! default value
999	call getin_p("alb_tho_eq",alb_tho_eq)
1000	if (is_master)write(,)trim(rname)//&
1001	" alb_tho_eq = ",alb_tho_eq
1002
1003	if (is_master)write(,)trim(rname)//&
1004	"emis tholins "
1005	emis_tho=1. ! default value
1006	call getin_p("emis_tho",emis_tho)
1007	if (is_master)write(,)trim(rname)//&
1008	" emis_tho = ",emis_tho
1009
1010	if (is_master)write(,)trim(rname)//&
1011	"emis tholins eq"
1012	emis_tho_eq=1. ! default value
1013	call getin_p("emis_tho_eq",emis_tho_eq)
1014	if (is_master)write(,)trim(rname)//&
1015	" emis_tho_eq = ",emis_tho_eq
1016
1017	if (is_master)write(,)trim(rname)//&
1018	"Latitude for diff albedo tholins"
1019	tholateq=25. ! default value
1020	call getin_p("tholateq",tholateq)
1021	if (is_master)write(,)trim(rname)//&
1022	" tholateq = ",tholateq
1023	tholatn=-1.
1024	tholats=0.
1025	tholone=0.
1026	tholonw=-1.
1027	alb_tho_spe=0.1 ! default value
1028	emis_tho_spe=1. ! default value
1029	call getin_p(" tholatn",tholatn)
1030	call getin_p(" tholats",tholats)
1031	call getin_p(" tholone",tholone)
1032	call getin_p(" tholonw",tholonw)
1033	if (is_master)write(,)trim(rname)//&
1034	" Values for special tholins albedo = ",tholatn,&
1035	tholats,tholone,tholonw,alb_tho_spe,emis_tho_spe
1036
1037	if (is_master)write(,)trim(rname)//&
1038	"Specific albedo"
1039	spelon1=-180. ! default value
1040	spelon2=180. ! default value
1041	spelat1=-90. ! default value
1042	spelat2=90. ! default value
1043	specalb=.false. ! default value
1044	if (is_master)write(,)trim(rname)//&
1045	"albedo/emis spe "
1046	albspe=0.1 ! default value
1047	emispe=1. ! default value
1048	call getin_p("spelon1",spelon1)
1049	call getin_p("spelon2",spelon2)
1050	call getin_p("spelat1",spelat1)
1051	call getin_p("spelat2",spelat2)
1052	call getin_p("specalb",specalb)
1053	call getin_p("albspe",albspe)
1054	call getin_p("emispe",emispe)
1055
1056	if (is_master)write(,)trim(rname)//&
1057	" specific = ",specalb
1058
1059	!******************** TI ***************************
1060
1061	if (is_master)write(,)trim(rname)//&
1062	"Change TI with time"
1063	changeti=.false. ! default value
1064	call getin_p("changeti",changeti)
1065	if (is_master)write(,)trim(rname)//&
1066	" changeti = ",changeti
1067	changetid=.false. ! default value for diurnal TI
1068	call getin_p("changetid",changetid)
1069	if (is_master)write(,)trim(rname)//&
1070	" changetid = ",changetid
1071
1072	if (is_master)write(,)trim(rname)//&
1073	"IT N2 "
1074	ITN2=800. ! default value
1075	call getin_p("ITN2",ITN2)
1076	if (is_master)write(,)trim(rname)//&
1077	" ITN2 = ",ITN2
1078	ITN2d=20. ! default value
1079	call getin_p("ITN2d",ITN2d)
1080	if (is_master)write(,)trim(rname)//&
1081	" ITN2d = ",ITN2d
1082
1083	if (is_master)write(,)trim(rname)//&
1084	"IT CH4"
1085	ITCH4=800. ! default value
1086	call getin_p("ITCH4",ITCH4)
1087	if (is_master)write(,)trim(rname)//&
1088	" ITCH4 = ",ITCH4
1089	ITCH4d=20. ! default value
1090	call getin_p("ITCH4d",ITCH4d)
1091	if (is_master)write(,)trim(rname)//&
1092	" ITCH4d = ",ITCH4d
1093
1094	if (is_master)write(,)trim(rname)//&
1095	"IT H2O"
1096	ITH2O=800. ! default value
1097	call getin_p("ITH2O",ITH2O)
1098	if (is_master)write(,)trim(rname)//&
1099	" ITH2O = ",ITH2O
1100	ITH2Od=20. ! default value
1101	call getin_p("ITH2Od",ITH2Od)
1102	if (is_master)write(,)trim(rname)//&
1103	" ITH2Od = ",ITH2Od
1104
1105	!************ COOLING *************
1106
1107	alpha_top=5e-11 ! default value
1108	call getin_p("alpha_top",alpha_top)
1109	if (is_master)write(,)trim(rname)//&
1110	" alpha_top = ",alpha_top
1111	pref=0.02 ! default value
1112	call getin_p("pref",pref)
1113	if (is_master)write(,)trim(rname)//&
1114	" pref = ",pref
1115	deltap=0.1 ! default value
1116	call getin_p("deltap",deltap)
1117	if (is_master)write(,)trim(rname)//&
1118	" deltap = ",deltap
1119
1120	!=================================
1121	! TWOLAY scheme and NH3 cloudare left for retrocompatibility only,
1122	! You should now use N-LAYER scheme (see below).
1123
1124	if (is_master) write(,)trim(rname)//&
1125	": Radiatively active two-layer aerosols?"
1126	aeroback2lay=.false. ! default value
1127	call getin_p("aeroback2lay",aeroback2lay)
1128	if (is_master) write(,)trim(rname)//": aeroback2lay = ",aeroback2lay
1129
1130	if (aeroback2lay.and.is_master) then
1131	print*,'Warning : The TWOLAY AEROSOL scheme is deprecated and buggy...'
1132	print*,'You should use the generic n-layer scheme (see aeronlay).'
1133	endif
1134
1135	if (is_master) write(,)trim(rname)//": Radiatively active ammonia cloud?"
1136	aeronh3=.false. ! default value
1137	call getin_p("aeronh3",aeronh3)
1138	if (is_master) write(,)trim(rname)//": aeronh3 = ",aeronh3
1139
1140	if (aeronh3.and.is_master) then
1141	print*,'Warning : You are using specific NH3 cloud scheme ...'
1142	print*,'You should use the generic n-layer scheme (see aeronlay).'
1143	endif
1144
1145	if (is_master) write(,)trim(rname)//&
1146	": TWOLAY AEROSOL: total optical depth "//&
1147	"in the tropospheric layer (visible)"
1148	obs_tau_col_tropo=8.D0
1149	call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
1150	if (is_master) write(,)trim(rname)//&
1151	": obs_tau_col_tropo = ",obs_tau_col_tropo
1152
1153	if (is_master) write(,)trim(rname)//&
1154	": TWOLAY AEROSOL: total optical depth "//&
1155	"in the stratospheric layer (visible)"
1156	obs_tau_col_strato=0.08D0
1157	call getin_p("obs_tau_col_strato",obs_tau_col_strato)
1158	if (is_master) write(,)trim(rname)//&
1159	": obs_tau_col_strato = ",obs_tau_col_strato
1160
1161	if (is_master) write(,)trim(rname)//&
1162	": TWOLAY AEROSOL: optprop_back2lay_vis?"
1163	optprop_back2lay_vis = 'optprop_saturn_vis_n20.dat'
1164	call getin_p("optprop_back2lay_vis",optprop_back2lay_vis)
1165	if (is_master) write(,)trim(rname)//&
1166	": optprop_back2lay_vis = ",trim(optprop_back2lay_vis)
1167
1168	if (is_master) write(,)trim(rname)//&
1169	": TWOLAY AEROSOL: optprop_back2lay_ir?"
1170	optprop_back2lay_ir = 'optprop_saturn_ir_n20.dat'
1171	call getin_p("optprop_back2lay_ir",optprop_back2lay_ir)
1172	if (is_master) write(,)trim(rname)//&
1173	": optprop_back2lay_ir = ",trim(optprop_back2lay_ir)
1174
1175	if (is_master) write(,)trim(rname)//&
1176	": TWOLAY AEROSOL: pres_bottom_tropo? in pa"
1177	pres_bottom_tropo=66000.0
1178	call getin_p("pres_bottom_tropo",pres_bottom_tropo)
1179	if (is_master) write(,)trim(rname)//&
1180	": pres_bottom_tropo = ",pres_bottom_tropo
1181
1182	if (is_master) write(,)trim(rname)//&
1183	": TWOLAY AEROSOL: pres_top_tropo? in pa"
1184	pres_top_tropo=18000.0
1185	call getin_p("pres_top_tropo",pres_top_tropo)
1186	if (is_master) write(,)trim(rname)//&
1187	": pres_top_tropo = ",pres_top_tropo
1188
1189	if (is_master) write(,)trim(rname)//&
1190	": TWOLAY AEROSOL: pres_bottom_strato? in pa"
1191	pres_bottom_strato=2000.0
1192	call getin_p("pres_bottom_strato",pres_bottom_strato)
1193	if (is_master) write(,)trim(rname)//&
1194	": pres_bottom_strato = ",pres_bottom_strato
1195
1196	! Sanity check
1197	if (pres_bottom_strato .gt. pres_top_tropo) then
1198	if(is_master) then
1199	print*,'Error : TWOLAY AEROSOL, Please ensure that in callphys.def'
1200	print*,'you have pres_top_tropo > pres_bottom_strato !'
1201	endif
1202	call abort_physic(rname," pres_top_tropo > pres_bottom_strato!",1)
1203	endif
1204
1205	if (is_master) write(,)trim(rname)//&
1206	": TWOLAY AEROSOL: pres_top_strato? in pa"
1207	pres_top_strato=100.0
1208	call getin_p("pres_top_strato",pres_top_strato)
1209	if (is_master) write(,)trim(rname)//&
1210	": pres_top_strato = ",pres_top_strato
1211
1212	if (is_master) write(,)trim(rname)//&
1213	": TWOLAY AEROSOL: particle size in the ", &
1214	"tropospheric layer, in meters"
1215	size_tropo=2.e-6
1216	call getin_p("size_tropo",size_tropo)
1217	if (is_master) write(,)trim(rname)//": size_tropo = ",size_tropo
1218
1219	if (is_master) write(,)trim(rname)//&
1220	": TWOLAY AEROSOL: particle size in the ", &
1221	"stratospheric layer, in meters"
1222	size_strato=1.e-7
1223	call getin_p("size_strato",size_strato)
1224	if (is_master) write(,)trim(rname)//": size_strato = ",size_strato
1225
1226	if (is_master) write(,)trim(rname)//&
1227	": NH3 (thin) cloud: total optical depth"
1228	tau_nh3_cloud=7.D0
1229	call getin_p("tau_nh3_cloud",tau_nh3_cloud)
1230	if (is_master) write(,)trim(rname)//": tau_nh3_cloud = ",tau_nh3_cloud
1231
1232	if (is_master) write(,)trim(rname)//": NH3 (thin) cloud pressure level"
1233	pres_nh3_cloud=7.D0
1234	call getin_p("pres_nh3_cloud",pres_nh3_cloud)
1235	if (is_master) write(,)trim(rname)//": pres_nh3_cloud = ",pres_nh3_cloud
1236
1237	if (is_master) write(,)trim(rname)//": NH3 (thin) cloud: particle sizes"
1238	size_nh3_cloud=3.e-6
1239	call getin_p("size_nh3_cloud",size_nh3_cloud)
1240	if (is_master) write(,)trim(rname)//": size_nh3_cloud = ",size_nh3_cloud
1241
1242	!=================================
1243	! Generic N-LAYER aerosol scheme
1244
1245	if (is_master) write(,)trim(rname)//&
1246	": Radiatively active generic n-layer aerosols?"
1247	aeronlay=.false. ! default value
1248	call getin_p("aeronlay",aeronlay)
1249	if (is_master) write(,)trim(rname)//": aeronlay = ",aeronlay
1250
1251	if (is_master) write(,)trim(rname)//&
1252	": Number of generic aerosols layers?"
1253	nlayaero=1 ! default value
1254	call getin_p("nlayaero",nlayaero)
1255	! Avoid to allocate arrays of size 0
1256	if (aeronlay .and. nlayaero.lt.1) then
1257	if (is_master) then
1258	print*, " You are trying to set no generic aerosols..."
1259	print*, " Set aeronlay=.false. instead ! I abort."
1260	endif
1261	call abort_physic(rname,"no generic aerosols...",1)
1262	endif
1263	if (.not. aeronlay) nlayaero=1
1264	if (is_master) write(,)trim(rname)//": nlayaero = ",nlayaero
1265
1266	! This is necessary, we just set the number of aerosol layers
1267	IF(.NOT.ALLOCATED(aeronlay_tauref)) ALLOCATE(aeronlay_tauref(nlayaero))
1268	IF(.NOT.ALLOCATED(aeronlay_lamref)) ALLOCATE(aeronlay_lamref(nlayaero))
1269	IF(.NOT.ALLOCATED(aeronlay_choice)) ALLOCATE(aeronlay_choice(nlayaero))
1270	IF(.NOT.ALLOCATED(aeronlay_pbot)) ALLOCATE(aeronlay_pbot(nlayaero))
1271	IF(.NOT.ALLOCATED(aeronlay_ptop)) ALLOCATE(aeronlay_ptop(nlayaero))
1272	IF(.NOT.ALLOCATED(aeronlay_sclhght)) ALLOCATE(aeronlay_sclhght(nlayaero))
1273	IF(.NOT.ALLOCATED(aeronlay_size)) ALLOCATE(aeronlay_size(nlayaero))
1274	IF(.NOT.ALLOCATED(aeronlay_nueff)) ALLOCATE(aeronlay_nueff(nlayaero))
1275	IF(.NOT.ALLOCATED(optprop_aeronlay_ir)) ALLOCATE(optprop_aeronlay_ir(nlayaero))
1276	IF(.NOT.ALLOCATED(optprop_aeronlay_vis)) ALLOCATE(optprop_aeronlay_vis(nlayaero))
1277
1278	if (is_master) write(,)trim(rname)//&
1279	": Generic n-layer aerosols: Optical depth at reference wavelenght"
1280	aeronlay_tauref=1.0E-1
1281	call getin_p("aeronlay_tauref",aeronlay_tauref)
1282	if (is_master) write(,)trim(rname)//&
1283	": aeronlay_tauref = ",aeronlay_tauref
1284
1285	if (is_master) write(,)trim(rname)//&
1286	": Generic n-layer aerosols: Reference wavelenght for optical depths (m)"
1287	aeronlay_lamref=0.6E-6
1288	call getin_p("aeronlay_lamref",aeronlay_lamref)
1289	if (is_master) write(,)trim(rname)//&
1290	": aeronlay_lamref = ",aeronlay_lamref
1291
1292	if (is_master) then
1293	write(,)trim(rname)//&
1294	": Generic n-layer aerosols: Vertical profile choice : "
1295	write(,)trim(rname)//&
1296	" (1) Tau btwn ptop and pbot follows atm. scale height"
1297	write(,)trim(rname)//&
1298	" or (2) Tau above pbot follows its own scale height"
1299	endif
1300	aeronlay_choice=1
1301	call getin_p("aeronlay_choice",aeronlay_choice)
1302	if (is_master) write(,)trim(rname)//&
1303	": aeronlay_choice = ",aeronlay_choice
1304
1305	if (is_master) write(,)trim(rname)//&
1306	": Generic n-layer aerosols: bottom pressures (Pa)"
1307	aeronlay_pbot=2000.0
1308	call getin_p("aeronlay_pbot",aeronlay_pbot)
1309	if (is_master) write(,)trim(rname)//": aeronlay_pbot = ",aeronlay_pbot
1310
1311	if (is_master) write(,)trim(rname)//&
1312	": Generic n-layer aerosols: (if choice=1) Top pressures (Pa) "
1313	aeronlay_ptop=300000.0
1314	call getin_p("aeronlay_ptop",aeronlay_ptop)
1315	if (is_master) write(,)trim(rname)//": aeronlay_ptop = ",aeronlay_ptop
1316
1317	if (is_master) write(,)trim(rname)//&
1318	": Generic n-layer aerosols: (if choice=2) Scale height / atm. scale height"
1319	aeronlay_sclhght=0.2
1320	call getin_p("aeronlay_sclhght",aeronlay_sclhght)
1321	if (is_master) write(,)trim(rname)//&
1322	": aeronlay_sclhght = ",aeronlay_sclhght
1323
1324	if (is_master) write(,)trim(rname)//&
1325	": Generic n-layer aerosols: particles effective radii (m)"
1326	aeronlay_size=1.e-6
1327	call getin_p("aeronlay_size",aeronlay_size)
1328	if (is_master) write(,)trim(rname)//": aeronlay_size = ",aeronlay_size
1329
1330	if (is_master) write(,)trim(rname)//&
1331	": Generic n-layer aerosols: particles radii effective variance"
1332	aeronlay_nueff=0.1
1333	call getin_p("aeronlay_nueff",aeronlay_nueff)
1334	if (is_master) write(,)trim(rname)//": aeronlay_nueff = ",aeronlay_nueff
1335
1336	if (is_master) write(,)trim(rname)//&
1337	": Generic n-layer aerosols: VIS optical properties file"
1338	optprop_aeronlay_vis = 'optprop_saturn_vis_n20.dat'
1339	call getin_p("optprop_aeronlay_vis",optprop_aeronlay_vis)
1340	if (is_master) write(,)trim(rname)//&
1341	": optprop_aeronlay_vis = ",optprop_aeronlay_vis
1342
1343	if (is_master) write(,)trim(rname)//&
1344	": Generic n-layer aerosols: IR optical properties file"
1345	optprop_aeronlay_ir = 'optprop_saturn_ir_n20.dat'
1346	call getin_p("optprop_aeronlay_ir",optprop_aeronlay_ir)
1347	if (is_master) write(,)trim(rname)//&
1348	": optprop_aeronlay_ir = ",optprop_aeronlay_ir
1349
1350
1351	!=================================
1352
1353	if (is_master) write(,)trim(rname)//&
1354	": Is the variable gas species radiatively active?"
1355	Tstrat=167.0
1356	varactive=.false.
1357	call getin_p("varactive",varactive)
1358	if (is_master) write(,)trim(rname)//": varactive = ",varactive
1359
1360	if (is_master) write(,)trim(rname)//&
1361	": Is the variable gas species distribution set?"
1362	varfixed=.false.
1363	call getin_p("varfixed",varfixed)
1364	if (is_master) write(,)trim(rname)//": varfixed = ",varfixed
1365
1366	if (is_master) write(,)trim(rname)//&
1367	": What is the saturation % of the variable species?"
1368	satval=0.8
1369	call getin_p("satval",satval)
1370	if (is_master) write(,)trim(rname)//": satval = ",satval
1371
1372
1373	! Test of incompatibility:
1374	! if varactive, then varfixed should be false
1375	if (varactive.and.varfixed) then
1376	call abort_physic(rname,'if varactive, varfixed must be OFF!',1)
1377	endif
1378
1379	if (is_master) write(,)trim(rname)//": Gravitationnal sedimentation ?"
1380	sedimentation=.false. ! default value
1381	call getin_p("sedimentation",sedimentation)
1382	if (is_master) write(,)trim(rname)//": sedimentation = ",sedimentation
1383
1384	if (is_master) write(,)trim(rname)//": Generic Condensation of tracers ?"
1385	generic_condensation=.false. !default value
1386	call getin_p("generic_condensation",generic_condensation)
1387	if (is_master) write(,)trim(rname)//": generic_condensation = ",generic_condensation
1388
1389	if (is_master) write(,)trim(rname)//": Spectral Dependant albedo ?"
1390	albedo_spectral_mode=.false. ! default value
1391	call getin_p("albedo_spectral_mode",albedo_spectral_mode)
1392	if (is_master) write(,)trim(rname)//&
1393	": albedo_spectral_mode = ",albedo_spectral_mode
1394
1395	if (is_master) then
1396	write(,)trim(rname)//": Snow albedo ?"
1397	write(,)trim(rname)//": If albedo_spectral_mode=.true., then this "
1398	write(,)trim(rname)//": corresponds to the 0.5 microns snow albedo."
1399	endif
1400	albedosnow=0.5 ! default value
1401	call getin_p("albedosnow",albedosnow)
1402	if (is_master) write(,)trim(rname)//": albedosnow = ",albedosnow
1403
1404	if (is_master) write(,)trim(rname)//": N2 ice albedo ?"
1405	albedon2ice=0.5 ! default value
1406	call getin_p("albedon2ice",albedon2ice)
1407	if (is_master) write(,)trim(rname)//": albedon2ice = ",albedon2ice
1408
1409	if (is_master) write(,)trim(rname)//": Maximum ice thickness ?"
1410	maxicethick=2.0 ! default value
1411	call getin_p("maxicethick",maxicethick)
1412	if (is_master) write(,)trim(rname)//": maxicethick = ",maxicethick
1413
1414	if (is_master) write(,)trim(rname)//": Minimum eddy mix coeff in 1D ?"
1415	kmixmin=1.0e-2 ! default value
1416	call getin_p("kmixmin",kmixmin)
1417	if (is_master) write(,)trim(rname)//": kmixmin = ",kmixmin
1418
1419	if (is_master) write(,)'Predefined Cp from dynamics is ',cpp,'J kg^-1 K^-1'
1420	if (is_master) write(,)'Predefined Mg from dynamics is ',mugaz,'amu'
1421
1422	force_cpp=.false. ! default value
1423	call getin_p("force_cpp",force_cpp)
1424	if (force_cpp) then
1425	if (is_master) write(,)trim(rname)//": force_cpp = ",force_cpp
1426	if (is_master) write(,)trim(rname)//": force_cpp is deprecated.",&
1427	"Set cpp_mugaz_mode=1 in callfis to emulate force_cpp=.true."
1428	call abort_physic(rname,"Anyway, you need to set force_cpp=.false. to continue.",1)
1429	endif
1430
1431	if (is_master) write(,)trim(rname)//&
1432	": where do you want your cpp/mugaz value to come from?",&
1433	"=> 0: dynamics (3d), 1: forced in callfis (1d), 2: computed from gases.def (1d)?"
1434	cpp_mugaz_mode = 0 ! default value
1435	call getin_p("cpp_mugaz_mode",cpp_mugaz_mode)
1436	if (is_master) write(,)trim(rname)//": cpp_mugaz_mode = ",cpp_mugaz_mode
1437
1438
1439	! Test of incompatibility:
1440
1441	if ((.not.tracer).and.(haze)) then
1442	call abort_physic(rname, 'if haze are on, tracers must be on!', 1)
1443
1444	endif
1445	if (callgasvis.and..not.callsoil) then
1446	call abort_physic(rname, 'if callgasvis is set, callsoil must be true', 1)
1447
1448	endif
1449	if (paleo.and..not.fast) then
1450	call abort_physic(rname, 'if paleo is set, fast must be true', 1)
1451
1452	endif
1453	if ((haze_proffix.or.haze_radproffix).and..not.aerohaze) then
1454	call abort_physic(rname, 'for now, haze/rad proffix only works w aerohaze=T', 1)
1455
1456	endif
1457	if (condcosurf.and.no_n2frost) then
1458	call abort_physic(rname, "CO surface condensation and no_n2frost are both active which may not be relevant", 1)
1459	end if
1460
1461	if ((cpp_mugaz_mode >= 1).and.(is_master).and.(ngrid>1)) then
1462	write(,)' !!!!! Be aware that having different values for cpp and mugaz in the dynamics and physics'
1463	write(,)' in 3D can result in very pathological behavior. You have been warned !!!!!'
1464	else if ((cpp_mugaz_mode == 0).and.(is_master).and.(ngrid==1)) then
1465	! for this specific 1D error we will remove run.def before aborting JL22
1466	call system("rm -rf run.def")
1467	call abort_physic(rname,"cpp_mugaz_mode must be >= 1 in 1d",1)
1468	endif
1469
1470	if (cpp_mugaz_mode == 1) then
1471	mugaz = -99999.
1472	if (is_master) write(,)trim(rname)//&
1473	": MEAN MOLECULAR MASS in g mol-1 ?"
1474	call getin_p("mugaz",mugaz)
1475	IF (mugaz.eq.-99999.) THEN
1476	call abort_physic(rname,"mugaz must be set if cpp_mugaz_mode = 1",1)
1477	ENDIF
1478	cpp = -99999.
1479	if (is_master) write(,)trim(rname)//&
1480	": SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
1481	call getin_p("cpp",cpp)
1482	IF (cpp.eq.-99999.) THEN
1483	PRINT *, "cpp must be set if cpp_mugaz_mode = 1"
1484	STOP
1485	ENDIF
1486	if (is_master) write(,)'New Cp from callfis is ',cpp,'J kg^-1 K^-1'
1487	if (is_master) write(,)'New Mg from callfis is ',mugaz,'amu'
1488
1489	endif ! of if (cpp_mugaz_mode == 1)
1490	call su_gases
1491	call calc_cpp_mugaz
1492
1493	if (is_master) then
1494	PRINT*,'--------------------------------------------'
1495	PRINT*
1496	PRINT*
1497	endif
1498	ELSE
1499	call abort_physic(rname,'Cannot read file callphys.def. Is it here ?',1)
1500	ENDIF ! of IF(iscallphys)
1501
1502	if (is_master) then
1503	PRINT*
1504	PRINT*,'inifis: daysec',daysec
1505	PRINT*
1506	PRINT*,'inifis: The radiative transfer is computed:'
1507	PRINT*,' each ',iradia,' physical time-step'
1508	PRINT,' or each ',iradiadtphys,' seconds'
1509	PRINT*
1510	endif
1511
1512	!-----------------------------------------------------------------------
1513	! Some more initialization:
1514	! ------------------------
1515
1516	! Initializations for comgeomfi_h
1517	#ifndef MESOSCALE
1518	totarea=SSUM(ngrid,parea,1)
1519	call planetwide_sumval(parea,totarea_planet)
1520
1521	!! those are defined in comdiurn_h.F90
1522	IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
1523	IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
1524	IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
1525	IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
1526
1527	DO ig=1,ngrid
1528	sinlat(ig)=sin(plat(ig))
1529	coslat(ig)=cos(plat(ig))
1530	sinlon(ig)=sin(plon(ig))
1531	coslon(ig)=cos(plon(ig))
1532	ENDDO
1533	#endif
1534	! initialize variables in radinc_h
1535	call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
1536
1537	! initialize variables and allocate arrays in radcommon_h
1538	call ini_radcommon_h(naerkind)
1539
1540	! allocate "comsoil_h" arrays
1541	call ini_comsoil_h(ngrid)
1542
1543	END SUBROUTINE inifis
1544
1545	END MODULE inifis_mod

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

Download in other formats:

Original Format