Context Navigation

inifis_mod.F90 @ 2803

Last change on this file since 2803 was 2803, checked in by aslmd, 2 years ago

Initialisation of Radiative Generic Condensable Aerosols

We can activate the scheme by putting aerogeneric = # of aerosols in callphys.def. This is the only needed thing for activating the radiative effects. They must be tracer in modern tracer.def

Commented out the abort if we use more than 4 aerosols

Added reading of optprop files for Radiative Generic Condensable Aerosols

We use the same file for IR and VI channel. For now, only MnS, Cr,Fe and Mg2SiO4 can be read. If you want to add another specie, check the code, it is explained how to quickly do that (right above the Initializations)

Added radii calculation for Radiative Generic Condensable aerosols

Changed the hardcoded size of the totalemit array

The hardcoded size is now 1900 instead of 100 so we don't exceed the array size when working at high spectral resolution (very rare case)

Added opacity computation for Radiative Generic Condensable aerosols

We do this computation in the same fashion as what's performed on water and dust.

switch iniaerosol and initracer order, to prepare for the RGCS scheme

Needed to switch the order of initialization so we can use the RGCS scheme without the assumption that ice and vap tracer of the same specie are following each other in the traceur.def file

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

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

Download in other formats:

Original Format