Context Navigation

inifis_mod.F90 @ 2937

Last change on this file since 2937 was 2831, checked in by emillour, 3 years ago
Generic PCM: Add the possibility to include Venus-like aerosols (triggered by option aerovenus=.true. in callphys.def); baseline is to use 5 distinct scatterers but each may be turned on/off (via aerovenus1, aerovenus2, aerovenus2p, aerovenus3, aerovenusUV flags which may be specified in callphys.def). GG
File size: 45.9 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(,)"Radiatively active Venus clouds ?"
690	aerovenus=.false. ! default value
691	call getin_p("aerovenus",aerovenus)
692	if (aerovenus) then
693	aerovenus1=.true. ! default value
694	aerovenus2=.true. ! default value
695	aerovenus2p=.true. ! default value
696	aerovenus3=.true. ! default value
697	aerovenusUV=.true. ! default value
698	else
699	aerovenus1=.false. ! default value
700	aerovenus2=.false. ! default value
701	aerovenus2p=.false. ! default value
702	aerovenus3=.false. ! default value
703	aerovenusUV=.false. ! default value
704	endif
705	! in case the user wants to specifically set/unset sub-options
706	call getin_p("aerovenus1",aerovenus1)
707	if (is_master) write(,)" aerovenus1 = ",aerovenus1
708	call getin_p("aerovenus2",aerovenus2)
709	if (is_master) write(,)" aerovenus2 = ",aerovenus2
710	call getin_p("aerovenus2p",aerovenus2p)
711	if (is_master) write(,)" aerovenus2p= ",aerovenus2p
712	call getin_p("aerovenus3",aerovenus3)
713	if (is_master) write(,)" aerovenus3 = ",aerovenus3
714	call getin_p("aerovenusUV",aerovenusUV)
715	if (is_master) write(,)" aerovenusUV= ",aerovenusUV
716	! Sanity check: if any of the aerovenus* is set to true
717	! then aeronovenus should also be set to true
718	if ((.not.aerovenus).and.(aerovenus1.or.aerovenus2.or.aerovenus2p.or.&
719	aerovenus3.or.aerovenusUV)) then
720	if(is_master) then
721	write(,)" Error, if you set some of the aerovenus* to true"
722	write(,)" then flag aerovenus should be set to true as well!"
723	endif
724	call abort_physic("inifis"," aerovenus* flags mismatch!",1)
725	endif
726
727	if (is_master) write(,)trim(rname)//&
728	": TWOLAY AEROSOL: total optical depth "//&
729	"in the tropospheric layer (visible)"
730	obs_tau_col_tropo=8.D0
731	call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
732	if (is_master) write(,)trim(rname)//&
733	": obs_tau_col_tropo = ",obs_tau_col_tropo
734
735	if (is_master) write(,)trim(rname)//&
736	": TWOLAY AEROSOL: total optical depth "//&
737	"in the stratospheric layer (visible)"
738	obs_tau_col_strato=0.08D0
739	call getin_p("obs_tau_col_strato",obs_tau_col_strato)
740	if (is_master) write(,)trim(rname)//&
741	": obs_tau_col_strato = ",obs_tau_col_strato
742
743	if (is_master) write(,)trim(rname)//&
744	": TWOLAY AEROSOL: optprop_back2lay_vis?"
745	optprop_back2lay_vis = 'optprop_saturn_vis_n20.dat'
746	call getin_p("optprop_back2lay_vis",optprop_back2lay_vis)
747	if (is_master) write(,)trim(rname)//&
748	": optprop_back2lay_vis = ",trim(optprop_back2lay_vis)
749
750	if (is_master) write(,)trim(rname)//&
751	": TWOLAY AEROSOL: optprop_back2lay_ir?"
752	optprop_back2lay_ir = 'optprop_saturn_ir_n20.dat'
753	call getin_p("optprop_back2lay_ir",optprop_back2lay_ir)
754	if (is_master) write(,)trim(rname)//&
755	": optprop_back2lay_ir = ",trim(optprop_back2lay_ir)
756
757	if (is_master) write(,)trim(rname)//&
758	": TWOLAY AEROSOL: pres_bottom_tropo? in pa"
759	pres_bottom_tropo=66000.0
760	call getin_p("pres_bottom_tropo",pres_bottom_tropo)
761	if (is_master) write(,)trim(rname)//&
762	": pres_bottom_tropo = ",pres_bottom_tropo
763
764	if (is_master) write(,)trim(rname)//&
765	": TWOLAY AEROSOL: pres_top_tropo? in pa"
766	pres_top_tropo=18000.0
767	call getin_p("pres_top_tropo",pres_top_tropo)
768	if (is_master) write(,)trim(rname)//&
769	": pres_top_tropo = ",pres_top_tropo
770
771	if (is_master) write(,)trim(rname)//&
772	": TWOLAY AEROSOL: pres_bottom_strato? in pa"
773	pres_bottom_strato=2000.0
774	call getin_p("pres_bottom_strato",pres_bottom_strato)
775	if (is_master) write(,)trim(rname)//&
776	": pres_bottom_strato = ",pres_bottom_strato
777
778	! Sanity check
779	if (pres_bottom_strato .gt. pres_top_tropo) then
780	if(is_master) then
781	print*,'Error : TWOLAY AEROSOL, Please ensure that in callphys.def'
782	print*,'you have pres_top_tropo > pres_bottom_strato !'
783	endif
784	call abort_physic(rname," pres_top_tropo > pres_bottom_strato!",1)
785	endif
786
787	if (is_master) write(,)trim(rname)//&
788	": TWOLAY AEROSOL: pres_top_strato? in pa"
789	pres_top_strato=100.0
790	call getin_p("pres_top_strato",pres_top_strato)
791	if (is_master) write(,)trim(rname)//&
792	": pres_top_strato = ",pres_top_strato
793
794	if (is_master) write(,)trim(rname)//&
795	": TWOLAY AEROSOL: particle size in the ", &
796	"tropospheric layer, in meters"
797	size_tropo=2.e-6
798	call getin_p("size_tropo",size_tropo)
799	if (is_master) write(,)trim(rname)//": size_tropo = ",size_tropo
800
801	if (is_master) write(,)trim(rname)//&
802	": TWOLAY AEROSOL: particle size in the ", &
803	"stratospheric layer, in meters"
804	size_strato=1.e-7
805	call getin_p("size_strato",size_strato)
806	if (is_master) write(,)trim(rname)//": size_strato = ",size_strato
807
808	if (is_master) write(,)trim(rname)//&
809	": NH3 (thin) cloud: total optical depth"
810	tau_nh3_cloud=7.D0
811	call getin_p("tau_nh3_cloud",tau_nh3_cloud)
812	if (is_master) write(,)trim(rname)//": tau_nh3_cloud = ",tau_nh3_cloud
813
814	if (is_master) write(,)trim(rname)//": NH3 (thin) cloud pressure level"
815	pres_nh3_cloud=7.D0
816	call getin_p("pres_nh3_cloud",pres_nh3_cloud)
817	if (is_master) write(,)trim(rname)//": pres_nh3_cloud = ",pres_nh3_cloud
818
819	if (is_master) write(,)trim(rname)//": NH3 (thin) cloud: particle sizes"
820	size_nh3_cloud=3.e-6
821	call getin_p("size_nh3_cloud",size_nh3_cloud)
822	if (is_master) write(,)trim(rname)//": size_nh3_cloud = ",size_nh3_cloud
823
824	!=================================
825	! Generic N-LAYER aerosol scheme
826
827	if (is_master) write(,)trim(rname)//&
828	": Radiatively active generic n-layer aerosols?"
829	aeronlay=.false. ! default value
830	call getin_p("aeronlay",aeronlay)
831	if (is_master) write(,)trim(rname)//": aeronlay = ",aeronlay
832
833	if (is_master) write(,)trim(rname)//&
834	": Number of generic aerosols layers?"
835	nlayaero=1 ! default value
836	call getin_p("nlayaero",nlayaero)
837	! Avoid to allocate arrays of size 0
838	if (aeronlay .and. nlayaero.lt.1) then
839	if (is_master) then
840	print*, " You are trying to set no generic aerosols..."
841	print*, " Set aeronlay=.false. instead ! I abort."
842	endif
843	call abort_physic(rname,"no generic aerosols...",1)
844	endif
845	if (.not. aeronlay) nlayaero=1
846	if (is_master) write(,)trim(rname)//": nlayaero = ",nlayaero
847
848	! This is necessary, we just set the number of aerosol layers
849	IF(.NOT.ALLOCATED(aeronlay_tauref)) ALLOCATE(aeronlay_tauref(nlayaero))
850	IF(.NOT.ALLOCATED(aeronlay_lamref)) ALLOCATE(aeronlay_lamref(nlayaero))
851	IF(.NOT.ALLOCATED(aeronlay_choice)) ALLOCATE(aeronlay_choice(nlayaero))
852	IF(.NOT.ALLOCATED(aeronlay_pbot)) ALLOCATE(aeronlay_pbot(nlayaero))
853	IF(.NOT.ALLOCATED(aeronlay_ptop)) ALLOCATE(aeronlay_ptop(nlayaero))
854	IF(.NOT.ALLOCATED(aeronlay_sclhght)) ALLOCATE(aeronlay_sclhght(nlayaero))
855	IF(.NOT.ALLOCATED(aeronlay_size)) ALLOCATE(aeronlay_size(nlayaero))
856	IF(.NOT.ALLOCATED(aeronlay_nueff)) ALLOCATE(aeronlay_nueff(nlayaero))
857	IF(.NOT.ALLOCATED(optprop_aeronlay_ir)) ALLOCATE(optprop_aeronlay_ir(nlayaero))
858	IF(.NOT.ALLOCATED(optprop_aeronlay_vis)) ALLOCATE(optprop_aeronlay_vis(nlayaero))
859
860	if (is_master) write(,)trim(rname)//&
861	": Generic n-layer aerosols: Optical depth at reference wavelenght"
862	aeronlay_tauref=1.0E-1
863	call getin_p("aeronlay_tauref",aeronlay_tauref)
864	if (is_master) write(,)trim(rname)//&
865	": aeronlay_tauref = ",aeronlay_tauref
866
867	if (is_master) write(,)trim(rname)//&
868	": Generic n-layer aerosols: Reference wavelenght for optical depths (m)"
869	aeronlay_lamref=0.6E-6
870	call getin_p("aeronlay_lamref",aeronlay_lamref)
871	if (is_master) write(,)trim(rname)//&
872	": aeronlay_lamref = ",aeronlay_lamref
873
874	if (is_master) then
875	write(,)trim(rname)//&
876	": Generic n-layer aerosols: Vertical profile choice : "
877	write(,)trim(rname)//&
878	" (1) Tau btwn ptop and pbot follows atm. scale height"
879	write(,)trim(rname)//&
880	" or (2) Tau above pbot follows its own scale height"
881	endif
882	aeronlay_choice=1
883	call getin_p("aeronlay_choice",aeronlay_choice)
884	if (is_master) write(,)trim(rname)//&
885	": aeronlay_choice = ",aeronlay_choice
886
887	if (is_master) write(,)trim(rname)//&
888	": Generic n-layer aerosols: bottom pressures (Pa)"
889	aeronlay_pbot=2000.0
890	call getin_p("aeronlay_pbot",aeronlay_pbot)
891	if (is_master) write(,)trim(rname)//": aeronlay_pbot = ",aeronlay_pbot
892
893	if (is_master) write(,)trim(rname)//&
894	": Generic n-layer aerosols: (if choice=1) Top pressures (Pa) "
895	aeronlay_ptop=300000.0
896	call getin_p("aeronlay_ptop",aeronlay_ptop)
897	if (is_master) write(,)trim(rname)//": aeronlay_ptop = ",aeronlay_ptop
898
899	if (is_master) write(,)trim(rname)//&
900	": Generic n-layer aerosols: (if choice=2) Scale height / atm. scale height"
901	aeronlay_sclhght=0.2
902	call getin_p("aeronlay_sclhght",aeronlay_sclhght)
903	if (is_master) write(,)trim(rname)//&
904	": aeronlay_sclhght = ",aeronlay_sclhght
905
906	if (is_master) write(,)trim(rname)//&
907	": Generic n-layer aerosols: particles effective radii (m)"
908	aeronlay_size=1.e-6
909	call getin_p("aeronlay_size",aeronlay_size)
910	if (is_master) write(,)trim(rname)//": aeronlay_size = ",aeronlay_size
911
912	if (is_master) write(,)trim(rname)//&
913	": Generic n-layer aerosols: particles radii effective variance"
914	aeronlay_nueff=0.1
915	call getin_p("aeronlay_nueff",aeronlay_nueff)
916	if (is_master) write(,)trim(rname)//": aeronlay_nueff = ",aeronlay_nueff
917
918	if (is_master) write(,)trim(rname)//&
919	": Generic n-layer aerosols: VIS optical properties file"
920	optprop_aeronlay_vis = 'optprop_saturn_vis_n20.dat'
921	call getin_p("optprop_aeronlay_vis",optprop_aeronlay_vis)
922	if (is_master) write(,)trim(rname)//&
923	": optprop_aeronlay_vis = ",optprop_aeronlay_vis
924
925	if (is_master) write(,)trim(rname)//&
926	": Generic n-layer aerosols: IR optical properties file"
927	optprop_aeronlay_ir = 'optprop_saturn_ir_n20.dat'
928	call getin_p("optprop_aeronlay_ir",optprop_aeronlay_ir)
929	if (is_master) write(,)trim(rname)//&
930	": optprop_aeronlay_ir = ",optprop_aeronlay_ir
931
932
933	!=================================
934
935	if (is_master) write(,)trim(rname)//&
936	": Cloud pressure level (with kastprof only):"
937	cloudlvl=0. ! default value
938	call getin_p("cloudlvl",cloudlvl)
939	if (is_master) write(,)trim(rname)//": cloudlvl = ",cloudlvl
940
941	if (is_master) write(,)trim(rname)//&
942	": Is the variable gas species radiatively active?"
943	Tstrat=167.0
944	varactive=.false.
945	call getin_p("varactive",varactive)
946	if (is_master) write(,)trim(rname)//": varactive = ",varactive
947
948	if (is_master) write(,)trim(rname)//&
949	": Is the variable gas species distribution set?"
950	varfixed=.false.
951	call getin_p("varfixed",varfixed)
952	if (is_master) write(,)trim(rname)//": varfixed = ",varfixed
953
954	if (is_master) write(,)trim(rname)//&
955	": What is the saturation % of the variable species?"
956	satval=0.8
957	call getin_p("satval",satval)
958	if (is_master) write(,)trim(rname)//": satval = ",satval
959
960
961	! Test of incompatibility:
962	! if varactive, then varfixed should be false
963	if (varactive.and.varfixed) then
964	call abort_physic(rname,'if varactive, varfixed must be OFF!',1)
965	endif
966
967	if (is_master) write(,)trim(rname)//": Gravitationnal sedimentation ?"
968	sedimentation=.false. ! default value
969	call getin_p("sedimentation",sedimentation)
970	if (is_master) write(,)trim(rname)//": sedimentation = ",sedimentation
971
972	if (is_master) write(,)trim(rname)//": Generic Condensation of tracers ?"
973	generic_condensation=.false. !default value
974	call getin_p("generic_condensation",generic_condensation)
975	if (is_master) write(,)trim(rname)//": generic_condensation = ",generic_condensation
976
977	if (is_master) write(,)trim(rname)//": Generic Condensation of tracers ?"
978	generic_rain=.false. !default value
979	call getin_p("generic_rain",generic_rain)
980	if (is_master) write(,)trim(rname)//": generic_rain = ",generic_rain
981
982	if (is_master) write(,)trim(rname)//": Compute water cycle ?"
983	water=.false. ! default value
984	call getin_p("water",water)
985	if (is_master) write(,)trim(rname)//": water = ",water
986
987	! Test of incompatibility:
988	! if water is true, there should be at least a tracer
989	if (water.and.(.not.tracer)) then
990	call abort_physic(rname,'if water is ON, tracer must be ON too!',1)
991	endif
992
993	if (is_master) write(,)trim(rname)//": Include water condensation ?"
994	watercond=.false. ! default value
995	call getin_p("watercond",watercond)
996	if (is_master) write(,)trim(rname)//": watercond = ",watercond
997
998	! Test of incompatibility:
999	! if watercond is used, then water should be used too
1000	if (watercond.and.(.not.water)) then
1001	call abort_physic(rname,&
1002	'if watercond is used, water should be used too',1)
1003	endif
1004
1005	if (is_master) write(,)trim(rname)//": Include water precipitation ?"
1006	waterrain=.false. ! default value
1007	call getin_p("waterrain",waterrain)
1008	if (is_master) write(,)trim(rname)//": waterrain = ",waterrain
1009
1010	if (is_master) write(,)trim(rname)//": Include surface hydrology ?"
1011	hydrology=.false. ! default value
1012	call getin_p("hydrology",hydrology)
1013	if (is_master) write(,)trim(rname)//": hydrology = ",hydrology
1014
1015	if (is_master) write(,)trim(rname)//": Evolve surface water sources ?"
1016	sourceevol=.false. ! default value
1017	call getin_p("sourceevol",sourceevol)
1018	if (is_master) write(,)trim(rname)//": sourceevol = ",sourceevol
1019
1020	if (is_master) write(,)trim(rname)//": Ice evolution timestep ?"
1021	icetstep=100.0 ! default value
1022	call getin_p("icetstep",icetstep)
1023	if (is_master) write(,)trim(rname)//": icetstep = ",icetstep
1024
1025	if (is_master) write(,)trim(rname)//": Spectral Dependant albedo ?"
1026	albedo_spectral_mode=.false. ! default value
1027	call getin_p("albedo_spectral_mode",albedo_spectral_mode)
1028	if (is_master) write(,)trim(rname)//&
1029	": albedo_spectral_mode = ",albedo_spectral_mode
1030
1031	if (is_master) then
1032	write(,)trim(rname)//": Snow albedo ?"
1033	write(,)trim(rname)//": If albedo_spectral_mode=.true., then this "
1034	write(,)trim(rname)//": corresponds to the 0.5 microns snow albedo."
1035	endif
1036	albedosnow=0.5 ! default value
1037	call getin_p("albedosnow",albedosnow)
1038	if (is_master) write(,)trim(rname)//": albedosnow = ",albedosnow
1039
1040	if (is_master) write(,)trim(rname)//": Ocean albedo ?"
1041	alb_ocean=0.07 ! default value
1042	call getin_p("alb_ocean",alb_ocean)
1043	if (is_master) write(,)trim(rname)//": alb_ocean = ",alb_ocean
1044
1045	if (is_master) write(,)trim(rname)//": CO2 ice albedo ?"
1046	albedoco2ice=0.5 ! default value
1047	call getin_p("albedoco2ice",albedoco2ice)
1048	if (is_master) write(,)trim(rname)//": albedoco2ice = ",albedoco2ice
1049
1050	if (is_master) write(,)trim(rname)//": Maximum ice thickness ?"
1051	maxicethick=2.0 ! default value
1052	call getin_p("maxicethick",maxicethick)
1053	if (is_master) write(,)trim(rname)//": maxicethick = ",maxicethick
1054
1055	if (is_master) write(,)trim(rname)//": Freezing point of seawater ?"
1056	Tsaldiff=-1.8 ! default value
1057	call getin_p("Tsaldiff",Tsaldiff)
1058	if (is_master) write(,)trim(rname)//": Tsaldiff = ",Tsaldiff
1059
1060	if (is_master) write(,)trim(rname)//": Minimum eddy mix coeff in 1D ?"
1061	kmixmin=1.0e-2 ! default value
1062	call getin_p("kmixmin",kmixmin)
1063	if (is_master) write(,)trim(rname)//": kmixmin = ",kmixmin
1064
1065	if (is_master) write(,)'Predefined Cp from dynamics is ',cpp,'J kg^-1 K^-1'
1066	if (is_master) write(,)'Predefined Mg from dynamics is ',mugaz,'amu'
1067
1068	force_cpp=.false. ! default value
1069	call getin_p("force_cpp",force_cpp)
1070	if (force_cpp) then
1071	if (is_master) write(,)trim(rname)//": force_cpp = ",force_cpp
1072	if (is_master) write(,)trim(rname)//": force_cpp is deprecated.",&
1073	"Set cpp_mugaz_mode=1 in callfis to emulate force_cpp=.true."
1074	call abort_physic(rname,"Anyway, you need to set force_cpp=.false. to continue.",1)
1075	endif
1076
1077	if (is_master) write(,)trim(rname)//&
1078	": where do you want your cpp/mugaz value to come from?",&
1079	"=> 0: dynamics (3d), 1: forced in callfis (1d), 2: computed from gases.def (1d)?"
1080	cpp_mugaz_mode = 0 ! default value
1081	call getin_p("cpp_mugaz_mode",cpp_mugaz_mode)
1082	if (is_master) write(,)trim(rname)//": cpp_mugaz_mode = ",cpp_mugaz_mode
1083
1084	if ((cpp_mugaz_mode >= 1).and.(is_master).and.(ngrid>1)) then
1085	write(,)' !!!!! Be aware that having different values for cpp and mugaz in the dynamics and physics'
1086	write(,)' in 3D can result in very pathological behavior. You have been warned !!!!!'
1087	else if ((cpp_mugaz_mode == 0).and.(is_master).and.(ngrid==1)) then
1088	! for this specific 1D error we will remove run.def before aborting JL22
1089	call system("rm -rf run.def")
1090	call abort_physic(rname,"cpp_mugaz_mode must be >= 1 in 1d",1)
1091	endif
1092
1093	if (cpp_mugaz_mode == 1) then
1094	mugaz = -99999.
1095	if (is_master) write(,)trim(rname)//&
1096	": MEAN MOLECULAR MASS in g mol-1 ?"
1097	call getin_p("mugaz",mugaz)
1098	IF (mugaz.eq.-99999.) THEN
1099	call abort_physic(rname,"mugaz must be set if cpp_mugaz_mode = 1",1)
1100	ENDIF
1101	cpp = -99999.
1102	if (is_master) write(,)trim(rname)//&
1103	": SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
1104	call getin_p("cpp",cpp)
1105	IF (cpp.eq.-99999.) THEN
1106	PRINT *, "cpp must be set if cpp_mugaz_mode = 1"
1107	STOP
1108	ENDIF
1109	if (is_master) write(,)'New Cp from callfis is ',cpp,'J kg^-1 K^-1'
1110	if (is_master) write(,)'New Mg from callfis is ',mugaz,'amu'
1111
1112	endif ! of if (cpp_mugaz_mode == 1)
1113	call su_gases
1114	call calc_cpp_mugaz
1115
1116	if (is_master) then
1117	PRINT*,'--------------------------------------------'
1118	PRINT*
1119	PRINT*
1120	endif
1121	ELSE
1122	call abort_physic(rname,'Cannot read file callphys.def. Is it here ?',1)
1123	ENDIF ! of IF(iscallphys)
1124
1125	if (is_master) then
1126	PRINT*
1127	PRINT*,'inifis: daysec',daysec
1128	PRINT*
1129	PRINT*,'inifis: The radiative transfer is computed:'
1130	PRINT*,' each ',iradia,' physical time-step'
1131	PRINT,' or each ',iradiadtphys,' seconds'
1132	PRINT*
1133	endif
1134
1135	!-----------------------------------------------------------------------
1136	! Some more initialization:
1137	! ------------------------
1138
1139	! Initializations for comgeomfi_h
1140	#ifndef MESOSCALE
1141	totarea=SSUM(ngrid,parea,1)
1142	call planetwide_sumval(parea,totarea_planet)
1143
1144	!! those are defined in comdiurn_h.F90
1145	IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
1146	IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
1147	IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
1148	IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
1149
1150	DO ig=1,ngrid
1151	sinlat(ig)=sin(plat(ig))
1152	coslat(ig)=cos(plat(ig))
1153	sinlon(ig)=sin(plon(ig))
1154	coslon(ig)=cos(plon(ig))
1155	ENDDO
1156	#endif
1157	! initialize variables in radinc_h
1158	call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
1159
1160	! allocate "comsoil_h" arrays
1161	call ini_comsoil_h(ngrid)
1162
1163	END SUBROUTINE inifis
1164
1165	END MODULE inifis_mod

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

Download in other formats:

Original Format