Context Navigation

inifis_mod.F90 @ 2664

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

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

Download in other formats:

Original Format