Context Navigation

inifis_mod.F90 @ 3046

Last change on this file since 3046 was 2972, checked in by emillour, 21 months ago
Generic PCM: Make number of scatterers fully dynamic (i.e. set in callphys.def and no longer by compilation option "-s #"). One should now specify naerkind = # in callphys.def (default is 0). EM
File size: 46.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_co2
15	use datafile_mod, only: datadir
16	use comdiurn_h, only: sinlat, coslat, sinlon, coslon
17	use comgeomfi_h, only: totarea, totarea_planet
18	use comsoil_h, only: ini_comsoil_h, nsoilmx, lay1_soil, alpha_soil
19	use time_phylmdz_mod, only: ecritphy,day_step,iphysiq, &
20	init_time, daysec, dtphys
21	use comcstfi_mod, only: rad, cpp, g, r, rcp, &
22	mugaz, pi, avocado
23	use planete_mod, only: nres
24	use planetwide_mod, only: planetwide_sumval
25	use callkeys_mod
26	use wstats_mod, only: callstats
27	use ioipsl_getin_p_mod, only : getin_p
28	use mod_phys_lmdz_para, only : is_parallel, is_master, bcast
29
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(,) trim(rname)//&
139	": Run with or without tracer transport ?"
140	tracer=.false. ! default value
141	call getin_p("tracer",tracer)
142	if (is_master) write(,) trim(rname)//": tracer = ",tracer
143
144	if (is_master) write(,) trim(rname)//&
145	": Run with or without atm mass update "//&
146	" due to tracer evaporation/condensation?"
147	mass_redistrib=.false. ! default value
148	call getin_p("mass_redistrib",mass_redistrib)
149	if (is_master) write(,) trim(rname)//&
150	": mass_redistrib = ",mass_redistrib
151
152	if (is_master) then
153	write(,) trim(rname)//": Diurnal cycle ?"
154	write(,) trim(rname)//&
155	": (if diurnal=false, diurnal averaged solar heating)"
156	endif
157	diurnal=.true. ! default value
158	call getin_p("diurnal",diurnal)
159	if (is_master) write(,) trim(rname)//": diurnal = ",diurnal
160
161	if (is_master) then
162	write(,) trim(rname)//": Seasonal cycle ?"
163	write(,) trim(rname)//&
164	": (if season=false, Ls stays constant, to value ", &
165	"set in 'start'"
166	endif
167	season=.true. ! default value
168	call getin_p("season",season)
169	if (is_master) write(,) trim(rname)//": season = ",season
170
171	if (is_master) write(,) trim(rname)//&
172	": No seasonal cycle: initial day to lock the run during restart"
173	noseason_day=0.0 ! default value
174	call getin_p("noseason_day",noseason_day)
175	if (is_master) write(,) trim(rname)//": noseason_day=", noseason_day
176
177	if (is_master) write(,) trim(rname)//": Tidally resonant rotation ?"
178	tlocked=.false. ! default value
179	call getin_p("tlocked",tlocked)
180	if (is_master) write(,) trim(rname)//": tlocked = ",tlocked
181
182	if (is_master) write(,) trim(rname)//": Saturn ring shadowing ?"
183	rings_shadow = .false.
184	call getin_p("rings_shadow", rings_shadow)
185	if (is_master) write(,) trim(rname)//": rings_shadow = ", rings_shadow
186
187	if (is_master) write(,) trim(rname)//&
188	": Compute latitude-dependent gravity field?"
189	oblate = .false.
190	call getin_p("oblate", oblate)
191	if (is_master) write(,) trim(rname)//": oblate = ", oblate
192
193	if (is_master) write(,) trim(rname)//": Flattening of the planet (a-b)/a "
194	flatten = 0.0
195	call getin_p("flatten", flatten)
196	if (is_master) write(,) trim(rname)//": flatten = ", flatten
197
198
199	if (is_master) write(,) trim(rname)//": Needed if oblate=.true.: J2"
200	J2 = 0.0
201	call getin_p("J2", J2)
202	if (is_master) write(,) trim(rname)//": J2 = ", J2
203
204	if (is_master) write(,) trim(rname)//&
205	": Needed if oblate=.true.: Planet mass (*1e24 kg)"
206	MassPlanet = 0.0
207	call getin_p("MassPlanet", MassPlanet)
208	if (is_master) write(,) trim(rname)//": MassPlanet = ", MassPlanet
209
210	if (is_master) write(,) trim(rname)//&
211	": Needed if oblate=.true.: Planet mean radius (m)"
212	Rmean = 0.0
213	call getin_p("Rmean", Rmean)
214	if (is_master) write(,) trim(rname)//": Rmean = ", Rmean
215
216	! Test of incompatibility:
217	! if tlocked, then diurnal should be false
218	if (tlocked.and.diurnal) then
219	call abort_physic(rname,"If diurnal=true, we should turn off tlocked",1)
220	endif
221
222	if (is_master) write(,) trim(rname)//": Tidal resonance ratio ?"
223	nres=0 ! default value
224	call getin_p("nres",nres)
225	if (is_master) write(,) trim(rname)//": nres = ",nres
226
227	if (is_master) write(,) trim(rname)//&
228	": Write some extra output to the screen ?"
229	lwrite=.false. ! default value
230	call getin_p("lwrite",lwrite)
231	if (is_master) write(,) trim(rname)//": lwrite = ",lwrite
232
233	if (is_master) write(,) trim(rname)//&
234	": Save statistics in file stats.nc ?"
235	callstats=.true. ! default value
236	call getin_p("callstats",callstats)
237	if (is_master) write(,) trim(rname)//": callstats = ",callstats
238
239	if (is_master) write(,) trim(rname)//&
240	": Test energy conservation of model physics ?"
241	enertest=.false. ! default value
242	call getin_p("enertest",enertest)
243	if (is_master) write(,) trim(rname)//": enertest = ",enertest
244
245	if (is_master) write(,) trim(rname)//": call radiative transfer ?"
246	callrad=.true. ! default value
247	call getin_p("callrad",callrad)
248	if (is_master) write(,) trim(rname)//": callrad = ",callrad
249
250	if (is_master) write(,) trim(rname)//&
251	": call correlated-k radiative transfer ?"
252	corrk=.true. ! default value
253	call getin_p("corrk",corrk)
254	if (is_master) write(,) trim(rname)//": corrk = ",corrk
255
256	if (is_master) write(,) trim(rname)//&
257	": prohibit calculations outside corrk T grid?"
258	strictboundcorrk=.true. ! default value
259	call getin_p("strictboundcorrk",strictboundcorrk)
260	if (is_master) write(,) trim(rname)//&
261	": strictboundcorrk = ",strictboundcorrk
262
263	if (is_master) write(,) trim(rname)//&
264	": prohibit calculations outside CIA T grid?"
265	strictboundcia=.true. ! default value
266	call getin_p("strictboundcia",strictboundcia)
267	if (is_master) write(,) trim(rname)//&
268	": strictboundcia = ",strictboundcia
269
270	if (is_master) write(,) trim(rname)//&
271	": Minimum atmospheric temperature for Planck function integration ?"
272	tplanckmin=30.0 ! default value
273	call getin_p("tplanckmin",tplanckmin)
274	if (is_master) write(,) trim(rname)//": tplanckmin = ",tplanckmin
275
276	if (is_master) write(,) trim(rname)//&
277	": Maximum atmospheric temperature for Planck function integration ?"
278	tplanckmax=1500.0 ! default value
279	call getin_p("tplanckmax",tplanckmax)
280	if (is_master) write(,) trim(rname)//": tplanckmax = ",tplanckmax
281
282	if (is_master) write(,) trim(rname)//&
283	": Temperature step for Planck function integration ?"
284	dtplanck=0.1 ! default value
285	call getin_p("dtplanck",dtplanck)
286	if (is_master) write(,) trim(rname)//": dtplanck = ",dtplanck
287
288	if (is_master) write(,) trim(rname)//&
289	": call gaseous absorption in the visible bands?"//&
290	" (matters only if callrad=T)"
291	callgasvis=.false. ! default value
292	call getin_p("callgasvis",callgasvis)
293	if (is_master) write(,) trim(rname)//": callgasvis = ",callgasvis
294
295	if (is_master) write(,) trim(rname)//&
296	": call continuum opacities in radiative transfer ?"//&
297	" (matters only if callrad=T)"
298	continuum=.true. ! default value
299	call getin_p("continuum",continuum)
300	if (is_master) write(,) trim(rname)//": continuum = ",continuum
301
302	if (is_master) write(,) trim(rname)//": version for H2H2 CIA file ?"
303	versH2H2cia=2011 ! default value (should be 2018 but retrocompatibility first)
304	call getin_p("versH2H2cia",versH2H2cia)
305	if (is_master) write(,) trim(rname)//": versH2H2cia = ",versH2H2cia
306	! Sanity check
307	if (versH2H2cia.ne.2011 .and. versH2H2cia.ne.2018) then
308	call abort_physic(rname,"Error: Choose a correct value (2011 or 2018) for versH2H2cia !",1)
309	endif
310
311	if (is_master) write(,) trim(rname)//&
312	": CIA - normal or equilibrium ortho-para mixture for H2?"
313	H2orthopara_mixture = 'normal'
314	call getin_p("H2orthopara_mixture",H2orthopara_mixture)
315	if (is_master) write(,)trim(rname)//&
316	": H2orthopara_mixture = ",trim(H2orthopara_mixture)
317
318	if (is_master) write(,) trim(rname)//&
319	": call turbulent vertical diffusion ?"
320	calldifv=.true. ! default value
321	call getin_p("calldifv",calldifv)
322	if (is_master) write(,) trim(rname)//": calldifv = ",calldifv
323
324	if (is_master) write(,) trim(rname)//": use turbdiff instead of vdifc ?"
325	UseTurbDiff=.true. ! default value
326	call getin_p("UseTurbDiff",UseTurbDiff)
327	if (is_master) write(,) trim(rname)//": UseTurbDiff = ",UseTurbDiff
328
329	if (is_master) write(,) trim(rname)//": call convective adjustment ?"
330	calladj=.true. ! default value
331	call getin_p("calladj",calladj)
332	if (is_master) write(,) trim(rname)//": calladj = ",calladj
333
334	if (is_master) write(,)trim(rname)//&
335	": call Lott's non-oro GWs parameterisation scheme ?"
336	calllott_nonoro=.false. ! default value
337	call getin_p("calllott_nonoro",calllott_nonoro)
338	if (is_master) write(,)trim(rname)//&
339	": calllott_nonoro = ",calllott_nonoro
340
341	if (is_master) write(,) trim(rname)//": call thermal plume model ?"
342	calltherm=.false. ! default value
343	call getin_p("calltherm",calltherm)
344	if (is_master) write(,) trim(rname)//": calltherm = ",calltherm
345
346	if (is_master) write(,) trim(rname)//": call CO2 condensation ?"
347	co2cond=.false. ! default value
348	call getin_p("co2cond",co2cond)
349	if (is_master) write(,) trim(rname)//": co2cond = ",co2cond
350	! Test of incompatibility
351	if (co2cond.and.(.not.tracer)) then
352	call abort_physic(rname,"Error: We need a CO2 ice tracer to condense CO2!",1)
353	endif
354
355	if (is_master) write(,) trim(rname)//": CO2 supersaturation level ?"
356	co2supsat=1.0 ! default value
357	call getin_p("co2supsat",co2supsat)
358	if (is_master) write(,) trim(rname)//": co2supsat = ",co2supsat
359
360	if (is_master) write(,) trim(rname)//&
361	": Radiative timescale for Newtonian cooling (in Earth days)?"
362	tau_relax=30. ! default value
363	call getin_p("tau_relax",tau_relax)
364	if (is_master) write(,) trim(rname)//": tau_relax = ",tau_relax
365	tau_relax=tau_relax243600 ! convert Earth days --> seconds
366
367	if (is_master) write(,)trim(rname)//&
368	": call thermal conduction in the soil ?"
369	callsoil=.true. ! default value
370	call getin_p("callsoil",callsoil)
371	if (is_master) write(,) trim(rname)//": callsoil = ",callsoil
372
373	if (is_master) write(,)trim(rname)//&
374	": Rad transfer is computed every iradia", &
375	" physical timestep"
376	iradia=1 ! default value
377	call getin_p("iradia",iradia)
378	if (is_master) write(,)trim(rname)//": iradia = ",iradia
379
380	if (is_master) write(,)trim(rname)//": Rayleigh scattering ?"
381	rayleigh=.false.
382	call getin_p("rayleigh",rayleigh)
383	if (is_master) write(,)trim(rname)//": rayleigh = ",rayleigh
384
385	if (is_master) write(,) trim(rname)//&
386	": Use blackbody for stellar spectrum ?"
387	stelbbody=.false. ! default value
388	call getin_p("stelbbody",stelbbody)
389	if (is_master) write(,) trim(rname)//": stelbbody = ",stelbbody
390
391	if (is_master) write(,) trim(rname)//": Stellar blackbody temperature ?"
392	stelTbb=5800.0 ! default value
393	call getin_p("stelTbb",stelTbb)
394	if (is_master) write(,) trim(rname)//": stelTbb = ",stelTbb
395
396	if (is_master) write(,) trim(rname)//": Output mean OLR in 1D?"
397	meanOLR=.false.
398	call getin_p("meanOLR",meanOLR)
399	if (is_master) write(,) trim(rname)//": meanOLR = ",meanOLR
400
401	if (is_master) write(,)trim(rname)//": Output spectral OLR in 3D?"
402	specOLR=.false.
403	call getin_p("specOLR",specOLR)
404	if (is_master) write(,)trim(rname)//": specOLR = ",specOLR
405
406	if (is_master) write(,)trim(rname)//&
407	": Output diagnostic optical thickness attenuation exp(-tau)?"
408	diagdtau=.false.
409	call getin_p("diagdtau",diagdtau)
410	if (is_master) write(,)trim(rname)//": diagdtau = ",diagdtau
411
412	if (is_master) write(,)trim(rname)//": Operate in kastprof mode?"
413	kastprof=.false.
414	call getin_p("kastprof",kastprof)
415	if (is_master) write(,)trim(rname)//": kastprof = ",kastprof
416
417	if (is_master) write(,)trim(rname)//&
418	": Uniform absorption in radiative transfer?"
419	graybody=.false.
420	call getin_p("graybody",graybody)
421	if (is_master) write(,)trim(rname)//": graybody = ",graybody
422
423	! Soil model
424	if (is_master) write(,)trim(rname)//&
425	": Number of sub-surface layers for soil scheme?"
426	! default value of nsoilmx set in comsoil_h
427	call getin_p("nsoilmx",nsoilmx)
428	if (is_master) write(,)trim(rname)//": nsoilmx=",nsoilmx
429
430	if (is_master) write(,)trim(rname)//&
431	": Thickness of topmost soil layer (m)?"
432	! default value of lay1_soil set in comsoil_h
433	call getin_p("lay1_soil",lay1_soil)
434	if (is_master) write(,)trim(rname)//": lay1_soil = ",lay1_soil
435
436	if (is_master) write(,)trim(rname)//&
437	": Coefficient for soil layer thickness distribution?"
438	! default value of alpha_soil set in comsoil_h
439	call getin_p("alpha_soil",alpha_soil)
440	if (is_master) write(,)trim(rname)//": alpha_soil = ",alpha_soil
441
442	! Slab Ocean
443	if (is_master) write(,) trim(rname)//": Use slab-ocean ?"
444	ok_slab_ocean=.false. ! default value
445	call getin_p("ok_slab_ocean",ok_slab_ocean)
446	if (is_master) write(,) trim(rname)//": ok_slab_ocean = ",ok_slab_ocean
447	! Sanity check: for now slab oncean only works in serial mode
448	if (ok_slab_ocean.and.is_parallel) then
449	call abort_physic(rname," Error: slab ocean should only be used in serial mode!",1)
450	endif
451
452	if (is_master) write(,) trim(rname)//": Use slab-sea-ice ?"
453	ok_slab_sic=.true. ! default value
454	call getin_p("ok_slab_sic",ok_slab_sic)
455	if (is_master) write(,) trim(rname)//": ok_slab_sic = ",ok_slab_sic
456
457	if (is_master) write(,) trim(rname)//&
458	": Use heat transport for the ocean ?"
459	ok_slab_heat_transp=.true. ! default value
460	call getin_p("ok_slab_heat_transp",ok_slab_heat_transp)
461	if (is_master) write(,) trim(rname)//&
462	": ok_slab_heat_transp = ",ok_slab_heat_transp
463
464	! Photochemistry and chemistry in the thermosphere
465
466	if (is_master) write(,) trim(rname)//": Use photochemistry ?"
467	photochem=.false. ! default value
468	call getin_p("photochem",photochem)
469	if (is_master) write(,) trim(rname)//": photochem = ",photochem
470
471	if (is_master) write(,) trim(rname)//": Use photolysis heat table ?"
472	photoheat=.false. ! default value
473	call getin_p("photoheat",photoheat)
474	if (is_master) write(,) trim(rname)//": photoheat = ",photoheat
475
476	if (is_master) write(,) trim(rname)//&
477	": Use photolysis online calculation ?"
478	jonline=.false. ! default value
479	call getin_p("jonline",jonline)
480	if (is_master) write(,) trim(rname)//": jonline = ",jonline
481
482	if (is_master) write(,) trim(rname)//": Use deposition ?"
483	depos=.false. ! default value
484	call getin_p("depos",depos)
485	if (is_master) write(,) trim(rname)//": depos = ",depos
486
487	if (is_master) write(,)trim(rname)//": Production of haze ?"
488	haze=.false. ! default value
489	call getin_p("haze",haze)
490	if (is_master) write(,)trim(rname)//": haze = ",haze
491
492	! Global1D mean and solar zenith angle
493
494	if (ngrid.eq.1) then
495	PRINT*, 'Simulate global averaged conditions ?'
496	global1d = .false. ! default value
497	call getin_p("global1d",global1d)
498	write(,) "global1d = ",global1d
499
500	! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
501	if (global1d.and.diurnal) then
502	call abort_physic(rname,'if global1d is true, diurnal must be set to false',1)
503	endif
504
505	if (global1d) then
506	PRINT *,'Solar Zenith angle (deg.) ?'
507	PRINT *,'(assumed for averaged solar flux S/4)'
508	szangle=60.0 ! default value
509	call getin_p("szangle",szangle)
510	write(,) "szangle = ",szangle
511	endif
512	endif ! of if (ngrid.eq.1)
513
514	! Test of incompatibility:
515	! if kastprof used, we must be in 1D
516	if (kastprof.and.(ngrid.gt.1)) then
517	call abort_physic(rname,'kastprof can only be used in 1D!',1)
518	endif
519
520	if (is_master) write(,)trim(rname)//&
521	": Stratospheric temperature for kastprof mode?"
522	Tstrat=167.0
523	call getin_p("Tstrat",Tstrat)
524	if (is_master) write(,)trim(rname)//": Tstrat = ",Tstrat
525
526	if (is_master) write(,)trim(rname)//": Remove lower boundary?"
527	nosurf=.false.
528	call getin_p("nosurf",nosurf)
529	if (is_master) write(,)trim(rname)//": nosurf = ",nosurf
530
531	! Tests of incompatibility:
532	if (nosurf.and.callsoil) then
533	if (is_master) then
534	write(,)trim(rname)//'nosurf not compatible with soil scheme!'
535	write(,)trim(rname)//'... got to make a choice!'
536	endif
537	call abort_physic(rname,"nosurf not compatible with soil scheme!",1)
538	endif
539
540	if (is_master) write(,)trim(rname)//": Add an internal heat flux?", &
541	"... matters only if callsoil=F"
542	intheat=0.
543	call getin_p("intheat",intheat)
544	if (is_master) write(,)trim(rname)//": intheat = ",intheat
545
546	if (is_master) write(,)trim(rname)//&
547	": Use Newtonian cooling for radiative transfer?"
548	newtonian=.false.
549	call getin_p("newtonian",newtonian)
550	if (is_master) write(,)trim(rname)//": newtonian = ",newtonian
551
552	! Tests of incompatibility:
553	if (newtonian.and.corrk) then
554	call abort_physic(rname,'newtonian not compatible with correlated-k!',1)
555	endif
556	if (newtonian.and.calladj) then
557	call abort_physic(rname,'newtonian not compatible with adjustment!',1)
558	endif
559	if (newtonian.and.calldifv) then
560	call abort_physic(rname,'newtonian not compatible with a boundary layer!',1)
561	endif
562
563	if (is_master) write(,)trim(rname)//": Test physics timescale in 1D?"
564	testradtimes=.false.
565	call getin_p("testradtimes",testradtimes)
566	if (is_master) write(,)trim(rname)//": testradtimes = ",testradtimes
567
568	! Test of incompatibility:
569	! if testradtimes used, we must be in 1D
570	if (testradtimes.and.(ngrid.gt.1)) then
571	call abort_physic(rname,'testradtimes can only be used in 1D!',1)
572	endif
573
574	if (is_master) write(,)trim(rname)//": Default planetary temperature?"
575	tplanet=215.0
576	call getin_p("tplanet",tplanet)
577	if (is_master) write(,)trim(rname)//": tplanet = ",tplanet
578
579	if (is_master) write(,)trim(rname)//": Which star?"
580	startype=1 ! default value = Sol
581	call getin_p("startype",startype)
582	if (is_master) write(,)trim(rname)//": startype = ",startype
583
584	if (is_master) write(,)trim(rname)//": Value of stellar flux at 1 AU?"
585	Fat1AU=1356.0 ! default value = Sol today
586	call getin_p("Fat1AU",Fat1AU)
587	if (is_master) write(,)trim(rname)//": Fat1AU = ",Fat1AU
588
589
590	! TRACERS:
591
592	if (is_master) write(,)trim(rname)//": Varying H2O cloud fraction?"
593	CLFvarying=.false. ! default value
594	call getin_p("CLFvarying",CLFvarying)
595	if (is_master) write(,)trim(rname)//": CLFvarying = ",CLFvarying
596
597	if (is_master) write(,)trim(rname)//&
598	": Value of fixed H2O cloud fraction?"
599	CLFfixval=1.0 ! default value
600	call getin_p("CLFfixval",CLFfixval)
601	if (is_master) write(,)trim(rname)//": CLFfixval = ",CLFfixval
602
603	if (is_master) write(,)trim(rname)//": fixed radii for Cloud particles?"
604	radfixed=.false. ! default value
605	call getin_p("radfixed",radfixed)
606	if (is_master) write(,)trim(rname)//": radfixed = ",radfixed
607
608	if(kastprof)then
609	radfixed=.true.
610	endif
611
612	if (is_master) write(,)trim(rname)//&
613	": Number mixing ratio of CO2 ice particles:"
614	Nmix_co2=1.e6 ! default value
615	call getin_p("Nmix_co2",Nmix_co2)
616	if (is_master) write(,)trim(rname)//": Nmix_co2 = ",Nmix_co2
617
618	if (is_master) write(,)trim(rname)//&
619	"Number of radiatively active aerosols:"
620	naerkind=0 ! default value
621	call getin_p("naerkind",naerkind)
622	if (is_master) write(,)trim(rname)//": naerkind = ",naerkind
623
624	if (is_master) write(,)trim(rname)//": Opacity of dust (if used):"
625	dusttau=0. ! default value
626	call getin_p("dusttau",dusttau)
627	if (is_master) write(,)trim(rname)//": dusttau = ",dusttau
628
629	if (is_master) write(,)trim(rname)//": Radiatively active CO2 aerosols?"
630	aeroco2=.false. ! default value
631	call getin_p("aeroco2",aeroco2)
632	if (is_master) write(,)trim(rname)//": aeroco2 = ",aeroco2
633
634	if (is_master) write(,)trim(rname)//": Fixed CO2 aerosol distribution?"
635	aerofixco2=.false. ! default value
636	call getin_p("aerofixco2",aerofixco2)
637	if (is_master) write(,)trim(rname)//": aerofixco2 = ",aerofixco2
638
639	if (is_master) write(,)trim(rname)//": Radiatively active water ice?"
640	aeroh2o=.false. ! default value
641	call getin_p("aeroh2o",aeroh2o)
642	if (is_master) write(,)trim(rname)//": aeroh2o = ",aeroh2o
643
644	if (is_master) write(,)trim(rname)//": Fixed H2O aerosol distribution?"
645	aerofixh2o=.false. ! default value
646	call getin_p("aerofixh2o",aerofixh2o)
647	if (is_master) write(,)trim(rname)//": aerofixh2o = ",aerofixh2o
648
649	if (is_master) write(,)trim(rname)//&
650	": Radiatively active sulfuric acid aerosols?"
651	aeroh2so4=.false. ! default value
652	call getin_p("aeroh2so4",aeroh2so4)
653	if (is_master) write(,)trim(rname)//": aeroh2so4 = ",aeroh2so4
654
655	if (is_master) write(,)trim(rname)//&
656	": Radiatively active auroral aerosols?"
657	aeroaurora=.false. ! default value
658	call getin_p("aeroaurora",aeroaurora)
659	if (is_master) write(,)trim(rname)//": aeroaurora = ",aeroaurora
660
661	if (is_master) write(,)trim(rname)//&
662	": Radiatively active Generic Condensable aerosols?"
663	aerogeneric=0 ! default value
664	call getin_p("aerogeneric",aerogeneric)
665	if (is_master) write(,)trim(rname)//":aerogeneric",aerogeneric
666
667
668	!=================================
669	! TWOLAY scheme and NH3 cloudare left for retrocompatibility only,
670	! You should now use N-LAYER scheme (see below).
671
672	if (is_master) write(,)trim(rname)//&
673	": Radiatively active two-layer aerosols?"
674	aeroback2lay=.false. ! default value
675	call getin_p("aeroback2lay",aeroback2lay)
676	if (is_master) write(,)trim(rname)//": aeroback2lay = ",aeroback2lay
677
678	if (aeroback2lay.and.is_master) then
679	print*,'Warning : The TWOLAY AEROSOL scheme is deprecated and buggy...'
680	print*,'You should use the generic n-layer scheme (see aeronlay).'
681	endif
682
683	if (is_master) write(,)trim(rname)//": Radiatively active ammonia cloud?"
684	aeronh3=.false. ! default value
685	call getin_p("aeronh3",aeronh3)
686	if (is_master) write(,)trim(rname)//": aeronh3 = ",aeronh3
687
688	if (aeronh3.and.is_master) then
689	print*,'Warning : You are using specific NH3 cloud scheme ...'
690	print*,'You should use the generic n-layer scheme (see aeronlay).'
691	endif
692
693	if (is_master) write(,)"Radiatively active Venus clouds ?"
694	aerovenus=.false. ! default value
695	call getin_p("aerovenus",aerovenus)
696	if (aerovenus) then
697	aerovenus1=.true. ! default value
698	aerovenus2=.true. ! default value
699	aerovenus2p=.true. ! default value
700	aerovenus3=.true. ! default value
701	aerovenusUV=.true. ! default value
702	else
703	aerovenus1=.false. ! default value
704	aerovenus2=.false. ! default value
705	aerovenus2p=.false. ! default value
706	aerovenus3=.false. ! default value
707	aerovenusUV=.false. ! default value
708	endif
709	! in case the user wants to specifically set/unset sub-options
710	call getin_p("aerovenus1",aerovenus1)
711	if (is_master) write(,)" aerovenus1 = ",aerovenus1
712	call getin_p("aerovenus2",aerovenus2)
713	if (is_master) write(,)" aerovenus2 = ",aerovenus2
714	call getin_p("aerovenus2p",aerovenus2p)
715	if (is_master) write(,)" aerovenus2p= ",aerovenus2p
716	call getin_p("aerovenus3",aerovenus3)
717	if (is_master) write(,)" aerovenus3 = ",aerovenus3
718	call getin_p("aerovenusUV",aerovenusUV)
719	if (is_master) write(,)" aerovenusUV= ",aerovenusUV
720	! Sanity check: if any of the aerovenus* is set to true
721	! then aeronovenus should also be set to true
722	if ((.not.aerovenus).and.(aerovenus1.or.aerovenus2.or.aerovenus2p.or.&
723	aerovenus3.or.aerovenusUV)) then
724	if(is_master) then
725	write(,)" Error, if you set some of the aerovenus* to true"
726	write(,)" then flag aerovenus should be set to true as well!"
727	endif
728	call abort_physic("inifis"," aerovenus* flags mismatch!",1)
729	endif
730
731	if (is_master) write(,)trim(rname)//&
732	": TWOLAY AEROSOL: total optical depth "//&
733	"in the tropospheric layer (visible)"
734	obs_tau_col_tropo=8.D0
735	call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
736	if (is_master) write(,)trim(rname)//&
737	": obs_tau_col_tropo = ",obs_tau_col_tropo
738
739	if (is_master) write(,)trim(rname)//&
740	": TWOLAY AEROSOL: total optical depth "//&
741	"in the stratospheric layer (visible)"
742	obs_tau_col_strato=0.08D0
743	call getin_p("obs_tau_col_strato",obs_tau_col_strato)
744	if (is_master) write(,)trim(rname)//&
745	": obs_tau_col_strato = ",obs_tau_col_strato
746
747	if (is_master) write(,)trim(rname)//&
748	": TWOLAY AEROSOL: optprop_back2lay_vis?"
749	optprop_back2lay_vis = 'optprop_saturn_vis_n20.dat'
750	call getin_p("optprop_back2lay_vis",optprop_back2lay_vis)
751	if (is_master) write(,)trim(rname)//&
752	": optprop_back2lay_vis = ",trim(optprop_back2lay_vis)
753
754	if (is_master) write(,)trim(rname)//&
755	": TWOLAY AEROSOL: optprop_back2lay_ir?"
756	optprop_back2lay_ir = 'optprop_saturn_ir_n20.dat'
757	call getin_p("optprop_back2lay_ir",optprop_back2lay_ir)
758	if (is_master) write(,)trim(rname)//&
759	": optprop_back2lay_ir = ",trim(optprop_back2lay_ir)
760
761	if (is_master) write(,)trim(rname)//&
762	": TWOLAY AEROSOL: pres_bottom_tropo? in pa"
763	pres_bottom_tropo=66000.0
764	call getin_p("pres_bottom_tropo",pres_bottom_tropo)
765	if (is_master) write(,)trim(rname)//&
766	": pres_bottom_tropo = ",pres_bottom_tropo
767
768	if (is_master) write(,)trim(rname)//&
769	": TWOLAY AEROSOL: pres_top_tropo? in pa"
770	pres_top_tropo=18000.0
771	call getin_p("pres_top_tropo",pres_top_tropo)
772	if (is_master) write(,)trim(rname)//&
773	": pres_top_tropo = ",pres_top_tropo
774
775	if (is_master) write(,)trim(rname)//&
776	": TWOLAY AEROSOL: pres_bottom_strato? in pa"
777	pres_bottom_strato=2000.0
778	call getin_p("pres_bottom_strato",pres_bottom_strato)
779	if (is_master) write(,)trim(rname)//&
780	": pres_bottom_strato = ",pres_bottom_strato
781
782	! Sanity check
783	if (pres_bottom_strato .gt. pres_top_tropo) then
784	if(is_master) then
785	print*,'Error : TWOLAY AEROSOL, Please ensure that in callphys.def'
786	print*,'you have pres_top_tropo > pres_bottom_strato !'
787	endif
788	call abort_physic(rname," pres_top_tropo > pres_bottom_strato!",1)
789	endif
790
791	if (is_master) write(,)trim(rname)//&
792	": TWOLAY AEROSOL: pres_top_strato? in pa"
793	pres_top_strato=100.0
794	call getin_p("pres_top_strato",pres_top_strato)
795	if (is_master) write(,)trim(rname)//&
796	": pres_top_strato = ",pres_top_strato
797
798	if (is_master) write(,)trim(rname)//&
799	": TWOLAY AEROSOL: particle size in the ", &
800	"tropospheric layer, in meters"
801	size_tropo=2.e-6
802	call getin_p("size_tropo",size_tropo)
803	if (is_master) write(,)trim(rname)//": size_tropo = ",size_tropo
804
805	if (is_master) write(,)trim(rname)//&
806	": TWOLAY AEROSOL: particle size in the ", &
807	"stratospheric layer, in meters"
808	size_strato=1.e-7
809	call getin_p("size_strato",size_strato)
810	if (is_master) write(,)trim(rname)//": size_strato = ",size_strato
811
812	if (is_master) write(,)trim(rname)//&
813	": NH3 (thin) cloud: total optical depth"
814	tau_nh3_cloud=7.D0
815	call getin_p("tau_nh3_cloud",tau_nh3_cloud)
816	if (is_master) write(,)trim(rname)//": tau_nh3_cloud = ",tau_nh3_cloud
817
818	if (is_master) write(,)trim(rname)//": NH3 (thin) cloud pressure level"
819	pres_nh3_cloud=7.D0
820	call getin_p("pres_nh3_cloud",pres_nh3_cloud)
821	if (is_master) write(,)trim(rname)//": pres_nh3_cloud = ",pres_nh3_cloud
822
823	if (is_master) write(,)trim(rname)//": NH3 (thin) cloud: particle sizes"
824	size_nh3_cloud=3.e-6
825	call getin_p("size_nh3_cloud",size_nh3_cloud)
826	if (is_master) write(,)trim(rname)//": size_nh3_cloud = ",size_nh3_cloud
827
828	!=================================
829	! Generic N-LAYER aerosol scheme
830
831	if (is_master) write(,)trim(rname)//&
832	": Radiatively active generic n-layer aerosols?"
833	aeronlay=.false. ! default value
834	call getin_p("aeronlay",aeronlay)
835	if (is_master) write(,)trim(rname)//": aeronlay = ",aeronlay
836
837	if (is_master) write(,)trim(rname)//&
838	": Number of generic aerosols layers?"
839	nlayaero=1 ! default value
840	call getin_p("nlayaero",nlayaero)
841	! Avoid to allocate arrays of size 0
842	if (aeronlay .and. nlayaero.lt.1) then
843	if (is_master) then
844	print*, " You are trying to set no generic aerosols..."
845	print*, " Set aeronlay=.false. instead ! I abort."
846	endif
847	call abort_physic(rname,"no generic aerosols...",1)
848	endif
849	if (.not. aeronlay) nlayaero=1
850	if (is_master) write(,)trim(rname)//": nlayaero = ",nlayaero
851
852	! This is necessary, we just set the number of aerosol layers
853	IF(.NOT.ALLOCATED(aeronlay_tauref)) ALLOCATE(aeronlay_tauref(nlayaero))
854	IF(.NOT.ALLOCATED(aeronlay_lamref)) ALLOCATE(aeronlay_lamref(nlayaero))
855	IF(.NOT.ALLOCATED(aeronlay_choice)) ALLOCATE(aeronlay_choice(nlayaero))
856	IF(.NOT.ALLOCATED(aeronlay_pbot)) ALLOCATE(aeronlay_pbot(nlayaero))
857	IF(.NOT.ALLOCATED(aeronlay_ptop)) ALLOCATE(aeronlay_ptop(nlayaero))
858	IF(.NOT.ALLOCATED(aeronlay_sclhght)) ALLOCATE(aeronlay_sclhght(nlayaero))
859	IF(.NOT.ALLOCATED(aeronlay_size)) ALLOCATE(aeronlay_size(nlayaero))
860	IF(.NOT.ALLOCATED(aeronlay_nueff)) ALLOCATE(aeronlay_nueff(nlayaero))
861	IF(.NOT.ALLOCATED(optprop_aeronlay_ir)) ALLOCATE(optprop_aeronlay_ir(nlayaero))
862	IF(.NOT.ALLOCATED(optprop_aeronlay_vis)) ALLOCATE(optprop_aeronlay_vis(nlayaero))
863
864	if (is_master) write(,)trim(rname)//&
865	": Generic n-layer aerosols: Optical depth at reference wavelenght"
866	aeronlay_tauref=1.0E-1
867	call getin_p("aeronlay_tauref",aeronlay_tauref)
868	if (is_master) write(,)trim(rname)//&
869	": aeronlay_tauref = ",aeronlay_tauref
870
871	if (is_master) write(,)trim(rname)//&
872	": Generic n-layer aerosols: Reference wavelenght for optical depths (m)"
873	aeronlay_lamref=0.6E-6
874	call getin_p("aeronlay_lamref",aeronlay_lamref)
875	if (is_master) write(,)trim(rname)//&
876	": aeronlay_lamref = ",aeronlay_lamref
877
878	if (is_master) then
879	write(,)trim(rname)//&
880	": Generic n-layer aerosols: Vertical profile choice : "
881	write(,)trim(rname)//&
882	" (1) Tau btwn ptop and pbot follows atm. scale height"
883	write(,)trim(rname)//&
884	" or (2) Tau above pbot follows its own scale height"
885	endif
886	aeronlay_choice=1
887	call getin_p("aeronlay_choice",aeronlay_choice)
888	if (is_master) write(,)trim(rname)//&
889	": aeronlay_choice = ",aeronlay_choice
890
891	if (is_master) write(,)trim(rname)//&
892	": Generic n-layer aerosols: bottom pressures (Pa)"
893	aeronlay_pbot=2000.0
894	call getin_p("aeronlay_pbot",aeronlay_pbot)
895	if (is_master) write(,)trim(rname)//": aeronlay_pbot = ",aeronlay_pbot
896
897	if (is_master) write(,)trim(rname)//&
898	": Generic n-layer aerosols: (if choice=1) Top pressures (Pa) "
899	aeronlay_ptop=300000.0
900	call getin_p("aeronlay_ptop",aeronlay_ptop)
901	if (is_master) write(,)trim(rname)//": aeronlay_ptop = ",aeronlay_ptop
902
903	if (is_master) write(,)trim(rname)//&
904	": Generic n-layer aerosols: (if choice=2) Scale height / atm. scale height"
905	aeronlay_sclhght=0.2
906	call getin_p("aeronlay_sclhght",aeronlay_sclhght)
907	if (is_master) write(,)trim(rname)//&
908	": aeronlay_sclhght = ",aeronlay_sclhght
909
910	if (is_master) write(,)trim(rname)//&
911	": Generic n-layer aerosols: particles effective radii (m)"
912	aeronlay_size=1.e-6
913	call getin_p("aeronlay_size",aeronlay_size)
914	if (is_master) write(,)trim(rname)//": aeronlay_size = ",aeronlay_size
915
916	if (is_master) write(,)trim(rname)//&
917	": Generic n-layer aerosols: particles radii effective variance"
918	aeronlay_nueff=0.1
919	call getin_p("aeronlay_nueff",aeronlay_nueff)
920	if (is_master) write(,)trim(rname)//": aeronlay_nueff = ",aeronlay_nueff
921
922	if (is_master) write(,)trim(rname)//&
923	": Generic n-layer aerosols: VIS optical properties file"
924	optprop_aeronlay_vis = 'optprop_saturn_vis_n20.dat'
925	call getin_p("optprop_aeronlay_vis",optprop_aeronlay_vis)
926	if (is_master) write(,)trim(rname)//&
927	": optprop_aeronlay_vis = ",optprop_aeronlay_vis
928
929	if (is_master) write(,)trim(rname)//&
930	": Generic n-layer aerosols: IR optical properties file"
931	optprop_aeronlay_ir = 'optprop_saturn_ir_n20.dat'
932	call getin_p("optprop_aeronlay_ir",optprop_aeronlay_ir)
933	if (is_master) write(,)trim(rname)//&
934	": optprop_aeronlay_ir = ",optprop_aeronlay_ir
935
936
937	!=================================
938
939	if (is_master) write(,)trim(rname)//&
940	": Cloud pressure level (with kastprof only):"
941	cloudlvl=0. ! default value
942	call getin_p("cloudlvl",cloudlvl)
943	if (is_master) write(,)trim(rname)//": cloudlvl = ",cloudlvl
944
945	if (is_master) write(,)trim(rname)//&
946	": Is the variable gas species radiatively active?"
947	Tstrat=167.0
948	varactive=.false.
949	call getin_p("varactive",varactive)
950	if (is_master) write(,)trim(rname)//": varactive = ",varactive
951
952	if (is_master) write(,)trim(rname)//&
953	": Is the variable gas species distribution set?"
954	varfixed=.false.
955	call getin_p("varfixed",varfixed)
956	if (is_master) write(,)trim(rname)//": varfixed = ",varfixed
957
958	if (is_master) write(,)trim(rname)//&
959	": What is the saturation % of the variable species?"
960	satval=0.8
961	call getin_p("satval",satval)
962	if (is_master) write(,)trim(rname)//": satval = ",satval
963
964
965	! Test of incompatibility:
966	! if varactive, then varfixed should be false
967	if (varactive.and.varfixed) then
968	call abort_physic(rname,'if varactive, varfixed must be OFF!',1)
969	endif
970
971	if (is_master) write(,)trim(rname)//": Gravitationnal sedimentation ?"
972	sedimentation=.false. ! default value
973	call getin_p("sedimentation",sedimentation)
974	if (is_master) write(,)trim(rname)//": sedimentation = ",sedimentation
975
976	if (is_master) write(,)trim(rname)//": Generic Condensation of tracers ?"
977	generic_condensation=.false. !default value
978	call getin_p("generic_condensation",generic_condensation)
979	if (is_master) write(,)trim(rname)//": generic_condensation = ",generic_condensation
980
981	if (is_master) write(,)trim(rname)//": Generic Condensation of tracers ?"
982	generic_rain=.false. !default value
983	call getin_p("generic_rain",generic_rain)
984	if (is_master) write(,)trim(rname)//": generic_rain = ",generic_rain
985
986	if (is_master) write(,)trim(rname)//": Compute water cycle ?"
987	water=.false. ! default value
988	call getin_p("water",water)
989	if (is_master) write(,)trim(rname)//": water = ",water
990
991	! Test of incompatibility:
992	! if water is true, there should be at least a tracer
993	if (water.and.(.not.tracer)) then
994	call abort_physic(rname,'if water is ON, tracer must be ON too!',1)
995	endif
996
997	if (is_master) write(,)trim(rname)//": Include water condensation ?"
998	watercond=.false. ! default value
999	call getin_p("watercond",watercond)
1000	if (is_master) write(,)trim(rname)//": watercond = ",watercond
1001
1002	! Test of incompatibility:
1003	! if watercond is used, then water should be used too
1004	if (watercond.and.(.not.water)) then
1005	call abort_physic(rname,&
1006	'if watercond is used, water should be used too',1)
1007	endif
1008
1009	if (is_master) write(,)trim(rname)//": Include water precipitation ?"
1010	waterrain=.false. ! default value
1011	call getin_p("waterrain",waterrain)
1012	if (is_master) write(,)trim(rname)//": waterrain = ",waterrain
1013
1014	if (is_master) write(,)trim(rname)//": Include surface hydrology ?"
1015	hydrology=.false. ! default value
1016	call getin_p("hydrology",hydrology)
1017	if (is_master) write(,)trim(rname)//": hydrology = ",hydrology
1018
1019	if (is_master) write(,)trim(rname)//": Evolve surface water sources ?"
1020	sourceevol=.false. ! default value
1021	call getin_p("sourceevol",sourceevol)
1022	if (is_master) write(,)trim(rname)//": sourceevol = ",sourceevol
1023
1024	if (is_master) write(,)trim(rname)//": Ice evolution timestep ?"
1025	icetstep=100.0 ! default value
1026	call getin_p("icetstep",icetstep)
1027	if (is_master) write(,)trim(rname)//": icetstep = ",icetstep
1028
1029	if (is_master) write(,)trim(rname)//": Spectral Dependant albedo ?"
1030	albedo_spectral_mode=.false. ! default value
1031	call getin_p("albedo_spectral_mode",albedo_spectral_mode)
1032	if (is_master) write(,)trim(rname)//&
1033	": albedo_spectral_mode = ",albedo_spectral_mode
1034
1035	if (is_master) then
1036	write(,)trim(rname)//": Snow albedo ?"
1037	write(,)trim(rname)//": If albedo_spectral_mode=.true., then this "
1038	write(,)trim(rname)//": corresponds to the 0.5 microns snow albedo."
1039	endif
1040	albedosnow=0.5 ! default value
1041	call getin_p("albedosnow",albedosnow)
1042	if (is_master) write(,)trim(rname)//": albedosnow = ",albedosnow
1043
1044	if (is_master) write(,)trim(rname)//": Ocean albedo ?"
1045	alb_ocean=0.07 ! default value
1046	call getin_p("alb_ocean",alb_ocean)
1047	if (is_master) write(,)trim(rname)//": alb_ocean = ",alb_ocean
1048
1049	if (is_master) write(,)trim(rname)//": CO2 ice albedo ?"
1050	albedoco2ice=0.5 ! default value
1051	call getin_p("albedoco2ice",albedoco2ice)
1052	if (is_master) write(,)trim(rname)//": albedoco2ice = ",albedoco2ice
1053
1054	if (is_master) write(,)trim(rname)//": Maximum ice thickness ?"
1055	maxicethick=2.0 ! default value
1056	call getin_p("maxicethick",maxicethick)
1057	if (is_master) write(,)trim(rname)//": maxicethick = ",maxicethick
1058
1059	if (is_master) write(,)trim(rname)//": Freezing point of seawater ?"
1060	Tsaldiff=-1.8 ! default value
1061	call getin_p("Tsaldiff",Tsaldiff)
1062	if (is_master) write(,)trim(rname)//": Tsaldiff = ",Tsaldiff
1063
1064	if (is_master) write(,)trim(rname)//": Minimum eddy mix coeff in 1D ?"
1065	kmixmin=1.0e-2 ! default value
1066	call getin_p("kmixmin",kmixmin)
1067	if (is_master) write(,)trim(rname)//": kmixmin = ",kmixmin
1068
1069	if (is_master) write(,)'Predefined Cp from dynamics is ',cpp,'J kg^-1 K^-1'
1070	if (is_master) write(,)'Predefined Mg from dynamics is ',mugaz,'amu'
1071
1072	force_cpp=.false. ! default value
1073	call getin_p("force_cpp",force_cpp)
1074	if (force_cpp) then
1075	if (is_master) write(,)trim(rname)//": force_cpp = ",force_cpp
1076	if (is_master) write(,)trim(rname)//": force_cpp is deprecated.",&
1077	"Set cpp_mugaz_mode=1 in callfis to emulate force_cpp=.true."
1078	call abort_physic(rname,"Anyway, you need to set force_cpp=.false. to continue.",1)
1079	endif
1080
1081	if (is_master) write(,)trim(rname)//&
1082	": where do you want your cpp/mugaz value to come from?",&
1083	"=> 0: dynamics (3d), 1: forced in callfis (1d), 2: computed from gases.def (1d)?"
1084	cpp_mugaz_mode = 0 ! default value
1085	call getin_p("cpp_mugaz_mode",cpp_mugaz_mode)
1086	if (is_master) write(,)trim(rname)//": cpp_mugaz_mode = ",cpp_mugaz_mode
1087
1088	if ((cpp_mugaz_mode >= 1).and.(is_master).and.(ngrid>1)) then
1089	write(,)' !!!!! Be aware that having different values for cpp and mugaz in the dynamics and physics'
1090	write(,)' in 3D can result in very pathological behavior. You have been warned !!!!!'
1091	else if ((cpp_mugaz_mode == 0).and.(is_master).and.(ngrid==1)) then
1092	! for this specific 1D error we will remove run.def before aborting JL22
1093	call system("rm -rf run.def")
1094	call abort_physic(rname,"cpp_mugaz_mode must be >= 1 in 1d",1)
1095	endif
1096
1097	if (cpp_mugaz_mode == 1) then
1098	mugaz = -99999.
1099	if (is_master) write(,)trim(rname)//&
1100	": MEAN MOLECULAR MASS in g mol-1 ?"
1101	call getin_p("mugaz",mugaz)
1102	IF (mugaz.eq.-99999.) THEN
1103	call abort_physic(rname,"mugaz must be set if cpp_mugaz_mode = 1",1)
1104	ENDIF
1105	cpp = -99999.
1106	if (is_master) write(,)trim(rname)//&
1107	": SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
1108	call getin_p("cpp",cpp)
1109	IF (cpp.eq.-99999.) THEN
1110	PRINT *, "cpp must be set if cpp_mugaz_mode = 1"
1111	STOP
1112	ENDIF
1113	if (is_master) write(,)'New Cp from callfis is ',cpp,'J kg^-1 K^-1'
1114	if (is_master) write(,)'New Mg from callfis is ',mugaz,'amu'
1115
1116	endif ! of if (cpp_mugaz_mode == 1)
1117	call su_gases
1118	call calc_cpp_mugaz
1119
1120	if (is_master) then
1121	PRINT*,'--------------------------------------------'
1122	PRINT*
1123	PRINT*
1124	endif
1125	ELSE
1126	call abort_physic(rname,'Cannot read file callphys.def. Is it here ?',1)
1127	ENDIF ! of IF(iscallphys)
1128
1129	if (is_master) then
1130	PRINT*
1131	PRINT*,'inifis: daysec',daysec
1132	PRINT*
1133	PRINT*,'inifis: The radiative transfer is computed:'
1134	PRINT*,' each ',iradia,' physical time-step'
1135	PRINT,' or each ',iradiadtphys,' seconds'
1136	PRINT*
1137	endif
1138
1139	!-----------------------------------------------------------------------
1140	! Some more initialization:
1141	! ------------------------
1142
1143	! Initializations for comgeomfi_h
1144	#ifndef MESOSCALE
1145	totarea=SSUM(ngrid,parea,1)
1146	call planetwide_sumval(parea,totarea_planet)
1147
1148	!! those are defined in comdiurn_h.F90
1149	IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
1150	IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
1151	IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
1152	IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
1153
1154	DO ig=1,ngrid
1155	sinlat(ig)=sin(plat(ig))
1156	coslat(ig)=cos(plat(ig))
1157	sinlon(ig)=sin(plon(ig))
1158	coslon(ig)=cos(plon(ig))
1159	ENDDO
1160	#endif
1161	! initialize variables in radinc_h
1162	call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
1163
1164	! initialize variables and allocate arrays in radcommon_h
1165	call ini_radcommon_h(naerkind)
1166
1167	! allocate "comsoil_h" arrays
1168	call ini_comsoil_h(ngrid)
1169
1170	END SUBROUTINE inifis
1171
1172	END MODULE inifis_mod

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

Download in other formats:

Original Format