Context Navigation

inifis_mod.F90 @ 2590

Last change on this file since 2590 was 2583, checked in by emillour, 4 years ago
Generic GCM: Bug fix (concerns OpenMP only) in inifis. And some cleanup in the output messages. EM
File size: 42.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)//&
243	": Check to see if cpp values used match gases.def ?"
244	check_cpp_match=.true. ! default value
245	call getin_p("check_cpp_match",check_cpp_match)
246	if (is_master) write(,) trim(rname)//&
247	": check_cpp_match = ",check_cpp_match
248
249	if (is_master) write(,) trim(rname)//": call radiative transfer ?"
250	callrad=.true. ! default value
251	call getin_p("callrad",callrad)
252	if (is_master) write(,) trim(rname)//": callrad = ",callrad
253
254	if (is_master) write(,) trim(rname)//&
255	": call correlated-k radiative transfer ?"
256	corrk=.true. ! default value
257	call getin_p("corrk",corrk)
258	if (is_master) write(,) trim(rname)//": corrk = ",corrk
259
260	if (is_master) write(,) trim(rname)//&
261	": prohibit calculations outside corrk T grid?"
262	strictboundcorrk=.true. ! default value
263	call getin_p("strictboundcorrk",strictboundcorrk)
264	if (is_master) write(,) trim(rname)//&
265	": strictboundcorrk = ",strictboundcorrk
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	": call turbulent vertical diffusion ?"
310	calldifv=.true. ! default value
311	call getin_p("calldifv",calldifv)
312	if (is_master) write(,) trim(rname)//": calldifv = ",calldifv
313
314	if (is_master) write(,) trim(rname)//": use turbdiff instead of vdifc ?"
315	UseTurbDiff=.true. ! default value
316	call getin_p("UseTurbDiff",UseTurbDiff)
317	if (is_master) write(,) trim(rname)//": UseTurbDiff = ",UseTurbDiff
318
319	if (is_master) write(,) trim(rname)//": call convective adjustment ?"
320	calladj=.true. ! default value
321	call getin_p("calladj",calladj)
322	if (is_master) write(,) trim(rname)//": calladj = ",calladj
323
324	if (is_master) write(,)trim(rname)//&
325	": call Lott's non-oro GWs parameterisation scheme ?"
326	calllott_nonoro=.false. ! default value
327	call getin_p("calllott_nonoro",calllott_nonoro)
328	if (is_master) write(,)trim(rname)//&
329	": calllott_nonoro = ",calllott_nonoro
330
331	if (is_master) write(,)trim(rname)//&
332	": Max Value of EP flux at launching altitude"
333	epflux_max=0.0 ! default value
334	call getin_p("epflux_max",epflux_max)
335	if (is_master) write(,)trim(rname)//": epflux_max = ",epflux_max
336
337	if (is_master) write(,) trim(rname)//": call thermal plume model ?"
338	calltherm=.false. ! default value
339	call getin_p("calltherm",calltherm)
340	if (is_master) write(,) trim(rname)//": calltherm = ",calltherm
341
342	if (is_master) write(,) trim(rname)//": call CO2 condensation ?"
343	co2cond=.false. ! default value
344	call getin_p("co2cond",co2cond)
345	if (is_master) write(,) trim(rname)//": co2cond = ",co2cond
346	! Test of incompatibility
347	if (co2cond.and.(.not.tracer)) then
348	call abort_physic(rname,"Error: We need a CO2 ice tracer to condense CO2!",1)
349	endif
350
351	if (is_master) write(,) trim(rname)//": CO2 supersaturation level ?"
352	co2supsat=1.0 ! default value
353	call getin_p("co2supsat",co2supsat)
354	if (is_master) write(,) trim(rname)//": co2supsat = ",co2supsat
355
356	if (is_master) write(,) trim(rname)//&
357	": Radiative timescale for Newtonian cooling (in Earth days)?"
358	tau_relax=30. ! default value
359	call getin_p("tau_relax",tau_relax)
360	if (is_master) write(,) trim(rname)//": tau_relax = ",tau_relax
361	tau_relax=tau_relax243600 ! convert Earth days --> seconds
362
363	if (is_master) write(,)trim(rname)//&
364	": call thermal conduction in the soil ?"
365	callsoil=.true. ! default value
366	call getin_p("callsoil",callsoil)
367	if (is_master) write(,) trim(rname)//": callsoil = ",callsoil
368
369	if (is_master) write(,)trim(rname)//&
370	": Rad transfer is computed every iradia", &
371	" physical timestep"
372	iradia=1 ! default value
373	call getin_p("iradia",iradia)
374	if (is_master) write(,)trim(rname)//": iradia = ",iradia
375
376	if (is_master) write(,)trim(rname)//": Rayleigh scattering ?"
377	rayleigh=.false.
378	call getin_p("rayleigh",rayleigh)
379	if (is_master) write(,)trim(rname)//": rayleigh = ",rayleigh
380
381	if (is_master) write(,) trim(rname)//&
382	": Use blackbody for stellar spectrum ?"
383	stelbbody=.false. ! default value
384	call getin_p("stelbbody",stelbbody)
385	if (is_master) write(,) trim(rname)//": stelbbody = ",stelbbody
386
387	if (is_master) write(,) trim(rname)//": Stellar blackbody temperature ?"
388	stelTbb=5800.0 ! default value
389	call getin_p("stelTbb",stelTbb)
390	if (is_master) write(,) trim(rname)//": stelTbb = ",stelTbb
391
392	if (is_master) write(,) trim(rname)//": Output mean OLR in 1D?"
393	meanOLR=.false.
394	call getin_p("meanOLR",meanOLR)
395	if (is_master) write(,) trim(rname)//": meanOLR = ",meanOLR
396
397	if (is_master) write(,)trim(rname)//": Output spectral OLR in 3D?"
398	specOLR=.false.
399	call getin_p("specOLR",specOLR)
400	if (is_master) write(,)trim(rname)//": specOLR = ",specOLR
401
402	if (is_master) write(,)trim(rname)//&
403	": Output diagnostic optical thickness attenuation exp(-tau)?"
404	diagdtau=.false.
405	call getin_p("diagdtau",diagdtau)
406	if (is_master) write(,)trim(rname)//": diagdtau = ",diagdtau
407
408	if (is_master) write(,)trim(rname)//": Operate in kastprof mode?"
409	kastprof=.false.
410	call getin_p("kastprof",kastprof)
411	if (is_master) write(,)trim(rname)//": kastprof = ",kastprof
412
413	if (is_master) write(,)trim(rname)//&
414	": Uniform absorption in radiative transfer?"
415	graybody=.false.
416	call getin_p("graybody",graybody)
417	if (is_master) write(,)trim(rname)//": graybody = ",graybody
418
419	! Soil model
420	if (is_master) write(,)trim(rname)//&
421	": Number of sub-surface layers for soil scheme?"
422	! default value of nsoilmx set in comsoil_h
423	call getin_p("nsoilmx",nsoilmx)
424	if (is_master) write(,)trim(rname)//": nsoilmx=",nsoilmx
425
426	if (is_master) write(,)trim(rname)//&
427	": Thickness of topmost soil layer (m)?"
428	! default value of lay1_soil set in comsoil_h
429	call getin_p("lay1_soil",lay1_soil)
430	if (is_master) write(,)trim(rname)//": lay1_soil = ",lay1_soil
431
432	if (is_master) write(,)trim(rname)//&
433	": Coefficient for soil layer thickness distribution?"
434	! default value of alpha_soil set in comsoil_h
435	call getin_p("alpha_soil",alpha_soil)
436	if (is_master) write(,)trim(rname)//": alpha_soil = ",alpha_soil
437
438	! Slab Ocean
439	if (is_master) write(,) trim(rname)//": Use slab-ocean ?"
440	ok_slab_ocean=.false. ! default value
441	call getin_p("ok_slab_ocean",ok_slab_ocean)
442	if (is_master) write(,) trim(rname)//": ok_slab_ocean = ",ok_slab_ocean
443	! Sanity check: for now slab oncean only works in serial mode
444	if (ok_slab_ocean.and.is_parallel) then
445	call abort_physic(rname," Error: slab ocean should only be used in serial mode!",1)
446	endif
447
448	if (is_master) write(,) trim(rname)//": Use slab-sea-ice ?"
449	ok_slab_sic=.true. ! default value
450	call getin_p("ok_slab_sic",ok_slab_sic)
451	if (is_master) write(,) trim(rname)//": ok_slab_sic = ",ok_slab_sic
452
453	if (is_master) write(,) trim(rname)//&
454	": Use heat transport for the ocean ?"
455	ok_slab_heat_transp=.true. ! default value
456	call getin_p("ok_slab_heat_transp",ok_slab_heat_transp)
457	if (is_master) write(,) trim(rname)//&
458	": ok_slab_heat_transp = ",ok_slab_heat_transp
459
460	! Photochemistry and chemistry in the thermosphere
461
462	if (is_master) write(,) trim(rname)//": Use photochemistry ?"
463	photochem=.false. ! default value
464	call getin_p("photochem",photochem)
465	if (is_master) write(,) trim(rname)//": photochem = ",photochem
466
467	if (is_master) write(,) trim(rname)//": Use photolysis heat table ?"
468	photoheat=.false. ! default value
469	call getin_p("photoheat",photoheat)
470	if (is_master) write(,) trim(rname)//": photoheat = ",photoheat
471
472	if (is_master) write(,) trim(rname)//&
473	": Use photolysis online calculation ?"
474	jonline=.false. ! default value
475	call getin_p("jonline",jonline)
476	if (is_master) write(,) trim(rname)//": jonline = ",jonline
477
478	if (is_master) write(,) trim(rname)//": Use deposition ?"
479	depos=.false. ! default value
480	call getin_p("depos",depos)
481	if (is_master) write(,) trim(rname)//": depos = ",depos
482
483	if (is_master) write(,)trim(rname)//": Production of haze ?"
484	haze=.false. ! default value
485	call getin_p("haze",haze)
486	if (is_master) write(,)trim(rname)//": haze = ",haze
487
488	! Global1D mean and solar zenith angle
489
490	if (ngrid.eq.1) then
491	PRINT*, 'Simulate global averaged conditions ?'
492	global1d = .false. ! default value
493	call getin_p("global1d",global1d)
494	write(,) "global1d = ",global1d
495
496	! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
497	if (global1d.and.diurnal) then
498	call abort_physic(rname,'if global1d is true, diurnal must be set to false',1)
499	endif
500
501	if (global1d) then
502	PRINT *,'Solar Zenith angle (deg.) ?'
503	PRINT *,'(assumed for averaged solar flux S/4)'
504	szangle=60.0 ! default value
505	call getin_p("szangle",szangle)
506	write(,) "szangle = ",szangle
507	endif
508	endif ! of if (ngrid.eq.1)
509
510	! Test of incompatibility:
511	! if kastprof used, we must be in 1D
512	if (kastprof.and.(ngrid.gt.1)) then
513	call abort_physic(rname,'kastprof can only be used in 1D!',1)
514	endif
515
516	if (is_master) write(,)trim(rname)//&
517	": Stratospheric temperature for kastprof mode?"
518	Tstrat=167.0
519	call getin_p("Tstrat",Tstrat)
520	if (is_master) write(,)trim(rname)//": Tstrat = ",Tstrat
521
522	if (is_master) write(,)trim(rname)//": Remove lower boundary?"
523	nosurf=.false.
524	call getin_p("nosurf",nosurf)
525	if (is_master) write(,)trim(rname)//": nosurf = ",nosurf
526
527	! Tests of incompatibility:
528	if (nosurf.and.callsoil) then
529	if (is_master) then
530	write(,)trim(rname)//'nosurf not compatible with soil scheme!'
531	write(,)trim(rname)//'... got to make a choice!'
532	endif
533	call abort_physic(rname,"nosurf not compatible with soil scheme!",1)
534	endif
535
536	if (is_master) write(,)trim(rname)//": Add an internal heat flux?", &
537	"... matters only if callsoil=F"
538	intheat=0.
539	call getin_p("intheat",intheat)
540	if (is_master) write(,)trim(rname)//": intheat = ",intheat
541
542	if (is_master) write(,)trim(rname)//&
543	": Use Newtonian cooling for radiative transfer?"
544	newtonian=.false.
545	call getin_p("newtonian",newtonian)
546	if (is_master) write(,)trim(rname)//": newtonian = ",newtonian
547
548	! Tests of incompatibility:
549	if (newtonian.and.corrk) then
550	call abort_physic(rname,'newtonian not compatible with correlated-k!',1)
551	endif
552	if (newtonian.and.calladj) then
553	call abort_physic(rname,'newtonian not compatible with adjustment!',1)
554	endif
555	if (newtonian.and.calldifv) then
556	call abort_physic(rname,'newtonian not compatible with a boundary layer!',1)
557	endif
558
559	if (is_master) write(,)trim(rname)//": Test physics timescale in 1D?"
560	testradtimes=.false.
561	call getin_p("testradtimes",testradtimes)
562	if (is_master) write(,)trim(rname)//": testradtimes = ",testradtimes
563
564	! Test of incompatibility:
565	! if testradtimes used, we must be in 1D
566	if (testradtimes.and.(ngrid.gt.1)) then
567	call abort_physic(rname,'testradtimes can only be used in 1D!',1)
568	endif
569
570	if (is_master) write(,)trim(rname)//": Default planetary temperature?"
571	tplanet=215.0
572	call getin_p("tplanet",tplanet)
573	if (is_master) write(,)trim(rname)//": tplanet = ",tplanet
574
575	if (is_master) write(,)trim(rname)//": Which star?"
576	startype=1 ! default value = Sol
577	call getin_p("startype",startype)
578	if (is_master) write(,)trim(rname)//": startype = ",startype
579
580	if (is_master) write(,)trim(rname)//": Value of stellar flux at 1 AU?"
581	Fat1AU=1356.0 ! default value = Sol today
582	call getin_p("Fat1AU",Fat1AU)
583	if (is_master) write(,)trim(rname)//": Fat1AU = ",Fat1AU
584
585
586	! TRACERS:
587
588	if (is_master) write(,)trim(rname)//": Varying H2O cloud fraction?"
589	CLFvarying=.false. ! default value
590	call getin_p("CLFvarying",CLFvarying)
591	if (is_master) write(,)trim(rname)//": CLFvarying = ",CLFvarying
592
593	if (is_master) write(,)trim(rname)//&
594	": Value of fixed H2O cloud fraction?"
595	CLFfixval=1.0 ! default value
596	call getin_p("CLFfixval",CLFfixval)
597	if (is_master) write(,)trim(rname)//": CLFfixval = ",CLFfixval
598
599	if (is_master) write(,)trim(rname)//": fixed radii for Cloud particles?"
600	radfixed=.false. ! default value
601	call getin_p("radfixed",radfixed)
602	if (is_master) write(,)trim(rname)//": radfixed = ",radfixed
603
604	if(kastprof)then
605	radfixed=.true.
606	endif
607
608	if (is_master) write(,)trim(rname)//&
609	": Number mixing ratio of CO2 ice particles:"
610	Nmix_co2=1.e6 ! default value
611	call getin_p("Nmix_co2",Nmix_co2)
612	if (is_master) write(,)trim(rname)//": Nmix_co2 = ",Nmix_co2
613
614	! write(,)"Number of radiatively active aerosols:"
615	! naerkind=0. ! default value
616	! call getin_p("naerkind",naerkind)
617	! write(,)" naerkind = ",naerkind
618
619	if (is_master) write(,)trim(rname)//": Opacity of dust (if used):"
620	dusttau=0. ! default value
621	call getin_p("dusttau",dusttau)
622	if (is_master) write(,)trim(rname)//": dusttau = ",dusttau
623
624	if (is_master) write(,)trim(rname)//": Radiatively active CO2 aerosols?"
625	aeroco2=.false. ! default value
626	call getin_p("aeroco2",aeroco2)
627	if (is_master) write(,)trim(rname)//": aeroco2 = ",aeroco2
628
629	if (is_master) write(,)trim(rname)//": Fixed CO2 aerosol distribution?"
630	aerofixco2=.false. ! default value
631	call getin_p("aerofixco2",aerofixco2)
632	if (is_master) write(,)trim(rname)//": aerofixco2 = ",aerofixco2
633
634	if (is_master) write(,)trim(rname)//": Radiatively active water ice?"
635	aeroh2o=.false. ! default value
636	call getin_p("aeroh2o",aeroh2o)
637	if (is_master) write(,)trim(rname)//": aeroh2o = ",aeroh2o
638
639	if (is_master) write(,)trim(rname)//": Fixed H2O aerosol distribution?"
640	aerofixh2o=.false. ! default value
641	call getin_p("aerofixh2o",aerofixh2o)
642	if (is_master) write(,)trim(rname)//": aerofixh2o = ",aerofixh2o
643
644	if (is_master) write(,)trim(rname)//&
645	": Radiatively active sulfuric acid aerosols?"
646	aeroh2so4=.false. ! default value
647	call getin_p("aeroh2so4",aeroh2so4)
648	if (is_master) write(,)trim(rname)//": aeroh2so4 = ",aeroh2so4
649
650	if (is_master) write(,)trim(rname)//&
651	": Radiatively active auroral aerosols?"
652	aeroaurora=.false. ! default value
653	call getin_p("aeroaurora",aeroaurora)
654	if (is_master) write(,)trim(rname)//": aeroaurora = ",aeroaurora
655
656	!=================================
657	! TWOLAY scheme and NH3 cloudare left for retrocompatibility only,
658	! You should now use N-LAYER scheme (see below).
659
660	if (is_master) write(,)trim(rname)//&
661	": Radiatively active two-layer aerosols?"
662	aeroback2lay=.false. ! default value
663	call getin_p("aeroback2lay",aeroback2lay)
664	if (is_master) write(,)trim(rname)//": aeroback2lay = ",aeroback2lay
665
666	if (aeroback2lay.and.is_master) then
667	print*,'Warning : The TWOLAY AEROSOL scheme is deprecated and buggy...'
668	print*,'You should use the generic n-layer scheme (see aeronlay).'
669	endif
670
671	if (is_master) write(,)trim(rname)//": Radiatively active ammonia cloud?"
672	aeronh3=.false. ! default value
673	call getin_p("aeronh3",aeronh3)
674	if (is_master) write(,)trim(rname)//": aeronh3 = ",aeronh3
675
676	if (aeronh3.and.is_master) then
677	print*,'Warning : You are using specific NH3 cloud scheme ...'
678	print*,'You should use the generic n-layer scheme (see aeronlay).'
679	endif
680
681	if (is_master) write(,)trim(rname)//&
682	": TWOLAY AEROSOL: total optical depth "//&
683	"in the tropospheric layer (visible)"
684	obs_tau_col_tropo=8.D0
685	call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
686	if (is_master) write(,)trim(rname)//&
687	": obs_tau_col_tropo = ",obs_tau_col_tropo
688
689	if (is_master) write(,)trim(rname)//&
690	": TWOLAY AEROSOL: total optical depth "//&
691	"in the stratospheric layer (visible)"
692	obs_tau_col_strato=0.08D0
693	call getin_p("obs_tau_col_strato",obs_tau_col_strato)
694	if (is_master) write(,)trim(rname)//&
695	": obs_tau_col_strato = ",obs_tau_col_strato
696
697	if (is_master) write(,)trim(rname)//&
698	": TWOLAY AEROSOL: optprop_back2lay_vis?"
699	optprop_back2lay_vis = 'optprop_saturn_vis_n20.dat'
700	call getin_p("optprop_back2lay_vis",optprop_back2lay_vis)
701	if (is_master) write(,)trim(rname)//&
702	": optprop_back2lay_vis = ",trim(optprop_back2lay_vis)
703
704	if (is_master) write(,)trim(rname)//&
705	": TWOLAY AEROSOL: optprop_back2lay_ir?"
706	optprop_back2lay_ir = 'optprop_saturn_ir_n20.dat'
707	call getin_p("optprop_back2lay_ir",optprop_back2lay_ir)
708	if (is_master) write(,)trim(rname)//&
709	": optprop_back2lay_ir = ",trim(optprop_back2lay_ir)
710
711	if (is_master) write(,)trim(rname)//&
712	": TWOLAY AEROSOL: pres_bottom_tropo? in pa"
713	pres_bottom_tropo=66000.0
714	call getin_p("pres_bottom_tropo",pres_bottom_tropo)
715	if (is_master) write(,)trim(rname)//&
716	": pres_bottom_tropo = ",pres_bottom_tropo
717
718	if (is_master) write(,)trim(rname)//&
719	": TWOLAY AEROSOL: pres_top_tropo? in pa"
720	pres_top_tropo=18000.0
721	call getin_p("pres_top_tropo",pres_top_tropo)
722	if (is_master) write(,)trim(rname)//&
723	": pres_top_tropo = ",pres_top_tropo
724
725	if (is_master) write(,)trim(rname)//&
726	": TWOLAY AEROSOL: pres_bottom_strato? in pa"
727	pres_bottom_strato=2000.0
728	call getin_p("pres_bottom_strato",pres_bottom_strato)
729	if (is_master) write(,)trim(rname)//&
730	": pres_bottom_strato = ",pres_bottom_strato
731
732	! Sanity check
733	if (pres_bottom_strato .gt. pres_top_tropo) then
734	if(is_master) then
735	print*,'Error : TWOLAY AEROSOL, Please ensure that in callphys.def'
736	print*,'you have pres_top_tropo > pres_bottom_strato !'
737	endif
738	call abort_physic(rname," pres_top_tropo > pres_bottom_strato!",1)
739	endif
740
741	if (is_master) write(,)trim(rname)//&
742	": TWOLAY AEROSOL: pres_top_strato? in pa"
743	pres_top_strato=100.0
744	call getin_p("pres_top_strato",pres_top_strato)
745	if (is_master) write(,)trim(rname)//&
746	": pres_top_strato = ",pres_top_strato
747
748	if (is_master) write(,)trim(rname)//&
749	": TWOLAY AEROSOL: particle size in the ", &
750	"tropospheric layer, in meters"
751	size_tropo=2.e-6
752	call getin_p("size_tropo",size_tropo)
753	if (is_master) write(,)trim(rname)//": size_tropo = ",size_tropo
754
755	if (is_master) write(,)trim(rname)//&
756	": TWOLAY AEROSOL: particle size in the ", &
757	"stratospheric layer, in meters"
758	size_strato=1.e-7
759	call getin_p("size_strato",size_strato)
760	if (is_master) write(,)trim(rname)//": size_strato = ",size_strato
761
762	if (is_master) write(,)trim(rname)//&
763	": NH3 (thin) cloud: total optical depth"
764	tau_nh3_cloud=7.D0
765	call getin_p("tau_nh3_cloud",tau_nh3_cloud)
766	if (is_master) write(,)trim(rname)//": tau_nh3_cloud = ",tau_nh3_cloud
767
768	if (is_master) write(,)trim(rname)//": NH3 (thin) cloud pressure level"
769	pres_nh3_cloud=7.D0
770	call getin_p("pres_nh3_cloud",pres_nh3_cloud)
771	if (is_master) write(,)trim(rname)//": pres_nh3_cloud = ",pres_nh3_cloud
772
773	if (is_master) write(,)trim(rname)//": NH3 (thin) cloud: particle sizes"
774	size_nh3_cloud=3.e-6
775	call getin_p("size_nh3_cloud",size_nh3_cloud)
776	if (is_master) write(,)trim(rname)//": size_nh3_cloud = ",size_nh3_cloud
777
778	!=================================
779	! Generic N-LAYER aerosol scheme
780
781	if (is_master) write(,)trim(rname)//&
782	": Radiatively active generic n-layer aerosols?"
783	aeronlay=.false. ! default value
784	call getin_p("aeronlay",aeronlay)
785	if (is_master) write(,)trim(rname)//": aeronlay = ",aeronlay
786
787	if (is_master) write(,)trim(rname)//&
788	": Number of generic aerosols layers?"
789	nlayaero=1 ! default value
790	call getin_p("nlayaero",nlayaero)
791	! Avoid to allocate arrays of size 0
792	if (aeronlay .and. nlayaero.lt.1) then
793	if (is_master) then
794	print*, " You are trying to set no generic aerosols..."
795	print*, " Set aeronlay=.false. instead ! I abort."
796	endif
797	call abort_physic(rname,"no generic aerosols...",1)
798	endif
799	if (.not. aeronlay) nlayaero=1
800	if (is_master) write(,)trim(rname)//": nlayaero = ",nlayaero
801
802	! This is necessary, we just set the number of aerosol layers
803	IF(.NOT.ALLOCATED(aeronlay_tauref)) ALLOCATE(aeronlay_tauref(nlayaero))
804	IF(.NOT.ALLOCATED(aeronlay_lamref)) ALLOCATE(aeronlay_lamref(nlayaero))
805	IF(.NOT.ALLOCATED(aeronlay_choice)) ALLOCATE(aeronlay_choice(nlayaero))
806	IF(.NOT.ALLOCATED(aeronlay_pbot)) ALLOCATE(aeronlay_pbot(nlayaero))
807	IF(.NOT.ALLOCATED(aeronlay_ptop)) ALLOCATE(aeronlay_ptop(nlayaero))
808	IF(.NOT.ALLOCATED(aeronlay_sclhght)) ALLOCATE(aeronlay_sclhght(nlayaero))
809	IF(.NOT.ALLOCATED(aeronlay_size)) ALLOCATE(aeronlay_size(nlayaero))
810	IF(.NOT.ALLOCATED(aeronlay_nueff)) ALLOCATE(aeronlay_nueff(nlayaero))
811	IF(.NOT.ALLOCATED(optprop_aeronlay_ir)) ALLOCATE(optprop_aeronlay_ir(nlayaero))
812	IF(.NOT.ALLOCATED(optprop_aeronlay_vis)) ALLOCATE(optprop_aeronlay_vis(nlayaero))
813
814	if (is_master) write(,)trim(rname)//&
815	": Generic n-layer aerosols: Optical depth at reference wavelenght"
816	aeronlay_tauref=1.0E-1
817	call getin_p("aeronlay_tauref",aeronlay_tauref)
818	if (is_master) write(,)trim(rname)//&
819	": aeronlay_tauref = ",aeronlay_tauref
820
821	if (is_master) write(,)trim(rname)//&
822	": Generic n-layer aerosols: Reference wavelenght for optical depths (m)"
823	aeronlay_lamref=0.6E-6
824	call getin_p("aeronlay_lamref",aeronlay_lamref)
825	if (is_master) write(,)trim(rname)//&
826	": aeronlay_lamref = ",aeronlay_lamref
827
828	if (is_master) then
829	write(,)trim(rname)//&
830	": Generic n-layer aerosols: Vertical profile choice : "
831	write(,)trim(rname)//&
832	" (1) Tau btwn ptop and pbot follows atm. scale height"
833	write(,)trim(rname)//&
834	" or (2) Tau above pbot follows its own scale height"
835	endif
836	aeronlay_choice=1
837	call getin_p("aeronlay_choice",aeronlay_choice)
838	if (is_master) write(,)trim(rname)//&
839	": aeronlay_choice = ",aeronlay_choice
840
841	if (is_master) write(,)trim(rname)//&
842	": Generic n-layer aerosols: bottom pressures (Pa)"
843	aeronlay_pbot=2000.0
844	call getin_p("aeronlay_pbot",aeronlay_pbot)
845	if (is_master) write(,)trim(rname)//": aeronlay_pbot = ",aeronlay_pbot
846
847	if (is_master) write(,)trim(rname)//&
848	": Generic n-layer aerosols: (if choice=1) Top pressures (Pa) "
849	aeronlay_ptop=300000.0
850	call getin_p("aeronlay_ptop",aeronlay_ptop)
851	if (is_master) write(,)trim(rname)//": aeronlay_ptop = ",aeronlay_ptop
852
853	if (is_master) write(,)trim(rname)//&
854	": Generic n-layer aerosols: (if choice=2) Scale height / atm. scale height"
855	aeronlay_sclhght=0.2
856	call getin_p("aeronlay_sclhght",aeronlay_sclhght)
857	if (is_master) write(,)trim(rname)//&
858	": aeronlay_sclhght = ",aeronlay_sclhght
859
860	if (is_master) write(,)trim(rname)//&
861	": Generic n-layer aerosols: particles effective radii (m)"
862	aeronlay_size=1.e-6
863	call getin_p("aeronlay_size",aeronlay_size)
864	if (is_master) write(,)trim(rname)//": aeronlay_size = ",aeronlay_size
865
866	if (is_master) write(,)trim(rname)//&
867	": Generic n-layer aerosols: particles radii effective variance"
868	aeronlay_nueff=0.1
869	call getin_p("aeronlay_nueff",aeronlay_nueff)
870	if (is_master) write(,)trim(rname)//": aeronlay_nueff = ",aeronlay_nueff
871
872	if (is_master) write(,)trim(rname)//&
873	": Generic n-layer aerosols: VIS optical properties file"
874	optprop_aeronlay_vis = 'optprop_saturn_vis_n20.dat'
875	call getin_p("optprop_aeronlay_vis",optprop_aeronlay_vis)
876	if (is_master) write(,)trim(rname)//&
877	": optprop_aeronlay_vis = ",optprop_aeronlay_vis
878
879	if (is_master) write(,)trim(rname)//&
880	": Generic n-layer aerosols: IR optical properties file"
881	optprop_aeronlay_ir = 'optprop_saturn_ir_n20.dat'
882	call getin_p("optprop_aeronlay_ir",optprop_aeronlay_ir)
883	if (is_master) write(,)trim(rname)//&
884	": optprop_aeronlay_ir = ",optprop_aeronlay_ir
885
886
887	!=================================
888
889	if (is_master) write(,)trim(rname)//&
890	": Cloud pressure level (with kastprof only):"
891	cloudlvl=0. ! default value
892	call getin_p("cloudlvl",cloudlvl)
893	if (is_master) write(,)trim(rname)//": cloudlvl = ",cloudlvl
894
895	if (is_master) write(,)trim(rname)//&
896	": Is the variable gas species radiatively active?"
897	Tstrat=167.0
898	varactive=.false.
899	call getin_p("varactive",varactive)
900	if (is_master) write(,)trim(rname)//": varactive = ",varactive
901
902	if (is_master) write(,)trim(rname)//&
903	": Is the variable gas species distribution set?"
904	varfixed=.false.
905	call getin_p("varfixed",varfixed)
906	if (is_master) write(,)trim(rname)//": varfixed = ",varfixed
907
908	if (is_master) write(,)trim(rname)//&
909	": What is the saturation % of the variable species?"
910	satval=0.8
911	call getin_p("satval",satval)
912	if (is_master) write(,)trim(rname)//": satval = ",satval
913
914
915	! Test of incompatibility:
916	! if varactive, then varfixed should be false
917	if (varactive.and.varfixed) then
918	call abort_physic(rname,'if varactive, varfixed must be OFF!',1)
919	endif
920
921	if (is_master) write(,)trim(rname)//": Gravitationnal sedimentation ?"
922	sedimentation=.false. ! default value
923	call getin_p("sedimentation",sedimentation)
924	if (is_master) write(,)trim(rname)//": sedimentation = ",sedimentation
925
926	if (is_master) write(,)trim(rname)//": Compute water cycle ?"
927	water=.false. ! default value
928	call getin_p("water",water)
929	if (is_master) write(,)trim(rname)//": water = ",water
930
931	! Test of incompatibility:
932	! if water is true, there should be at least a tracer
933	if (water.and.(.not.tracer)) then
934	call abort_physic(rname,'if water is ON, tracer must be ON too!',1)
935	endif
936
937	if (is_master) write(,)trim(rname)//": Include water condensation ?"
938	watercond=.false. ! default value
939	call getin_p("watercond",watercond)
940	if (is_master) write(,)trim(rname)//": watercond = ",watercond
941
942	! Test of incompatibility:
943	! if watercond is used, then water should be used too
944	if (watercond.and.(.not.water)) then
945	call abort_physic(rname,&
946	'if watercond is used, water should be used too',1)
947	endif
948
949	if (is_master) write(,)trim(rname)//": Include water precipitation ?"
950	waterrain=.false. ! default value
951	call getin_p("waterrain",waterrain)
952	if (is_master) write(,)trim(rname)//": waterrain = ",waterrain
953
954	if (is_master) write(,)trim(rname)//": Include surface hydrology ?"
955	hydrology=.false. ! default value
956	call getin_p("hydrology",hydrology)
957	if (is_master) write(,)trim(rname)//": hydrology = ",hydrology
958
959	if (is_master) write(,)trim(rname)//": Evolve surface water sources ?"
960	sourceevol=.false. ! default value
961	call getin_p("sourceevol",sourceevol)
962	if (is_master) write(,)trim(rname)//": sourceevol = ",sourceevol
963
964	if (is_master) write(,)trim(rname)//": Ice evolution timestep ?"
965	icetstep=100.0 ! default value
966	call getin_p("icetstep",icetstep)
967	if (is_master) write(,)trim(rname)//": icetstep = ",icetstep
968
969	if (is_master) write(,)trim(rname)//": Spectral Dependant albedo ?"
970	albedo_spectral_mode=.false. ! default value
971	call getin_p("albedo_spectral_mode",albedo_spectral_mode)
972	if (is_master) write(,)trim(rname)//&
973	": albedo_spectral_mode = ",albedo_spectral_mode
974
975	if (is_master) then
976	write(,)trim(rname)//": Snow albedo ?"
977	write(,)trim(rname)//": If albedo_spectral_mode=.true., then this "
978	write(,)trim(rname)//": corresponds to the 0.5 microns snow albedo."
979	endif
980	albedosnow=0.5 ! default value
981	call getin_p("albedosnow",albedosnow)
982	if (is_master) write(,)trim(rname)//": albedosnow = ",albedosnow
983
984	if (is_master) write(,)trim(rname)//": Ocean albedo ?"
985	alb_ocean=0.07 ! default value
986	call getin_p("alb_ocean",alb_ocean)
987	if (is_master) write(,)trim(rname)//": alb_ocean = ",alb_ocean
988
989	if (is_master) write(,)trim(rname)//": CO2 ice albedo ?"
990	albedoco2ice=0.5 ! default value
991	call getin_p("albedoco2ice",albedoco2ice)
992	if (is_master) write(,)trim(rname)//": albedoco2ice = ",albedoco2ice
993
994	if (is_master) write(,)trim(rname)//": Maximum ice thickness ?"
995	maxicethick=2.0 ! default value
996	call getin_p("maxicethick",maxicethick)
997	if (is_master) write(,)trim(rname)//": maxicethick = ",maxicethick
998
999	if (is_master) write(,)trim(rname)//": Freezing point of seawater ?"
1000	Tsaldiff=-1.8 ! default value
1001	call getin_p("Tsaldiff",Tsaldiff)
1002	if (is_master) write(,)trim(rname)//": Tsaldiff = ",Tsaldiff
1003
1004	if (is_master) write(,)trim(rname)//": Minimum eddy mix coeff in 1D ?"
1005	kmixmin=1.0e-2 ! default value
1006	call getin_p("kmixmin",kmixmin)
1007	if (is_master) write(,)trim(rname)//": kmixmin = ",kmixmin
1008
1009	if (is_master) write(,)trim(rname)//&
1010	": Does user want to force cpp and mugaz?"
1011	force_cpp=.false. ! default value
1012	call getin_p("force_cpp",force_cpp)
1013	if (is_master) write(,)trim(rname)//": force_cpp = ",force_cpp
1014
1015	if (force_cpp) then
1016	mugaz = -99999.
1017	if (is_master) write(,)trim(rname)//&
1018	": MEAN MOLECULAR MASS in g mol-1 ?"
1019	call getin_p("mugaz",mugaz)
1020	IF (mugaz.eq.-99999.) THEN
1021	call abort_physic(rname,"mugaz must be set if force_cpp = T",1)
1022	ELSE
1023	if (is_master) write(,)trim(rname)//": mugaz=",mugaz
1024	ENDIF
1025	cpp = -99999.
1026	if (is_master) write(,)trim(rname)//&
1027	": SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?"
1028	call getin_p("cpp",cpp)
1029	IF (cpp.eq.-99999.) THEN
1030	PRINT *, "cpp must be set if force_cpp = T"
1031	STOP
1032	ELSE
1033	if (is_master) write(,)trim(rname)//": cpp=",cpp
1034	ENDIF
1035	endif ! of if (force_cpp)
1036	call su_gases
1037	call calc_cpp_mugaz
1038
1039	if (is_master) then
1040	PRINT*,'--------------------------------------------'
1041	PRINT*
1042	PRINT*
1043	endif
1044	ELSE
1045	call abort_physic(rname,'Cannot read file callphys.def. Is it here ?',1)
1046	ENDIF ! of IF(iscallphys)
1047
1048	if (is_master) then
1049	PRINT*
1050	PRINT*,'inifis: daysec',daysec
1051	PRINT*
1052	PRINT*,'inifis: The radiative transfer is computed:'
1053	PRINT*,' each ',iradia,' physical time-step'
1054	PRINT,' or each ',iradiadtphys,' seconds'
1055	PRINT*
1056	endif
1057
1058	!-----------------------------------------------------------------------
1059	! Some more initialization:
1060	! ------------------------
1061
1062	! Initializations for comgeomfi_h
1063	#ifndef MESOSCALE
1064	totarea=SSUM(ngrid,parea,1)
1065	call planetwide_sumval(parea,totarea_planet)
1066
1067	!! those are defined in comdiurn_h.F90
1068	IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
1069	IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
1070	IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
1071	IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
1072
1073	DO ig=1,ngrid
1074	sinlat(ig)=sin(plat(ig))
1075	coslat(ig)=cos(plat(ig))
1076	sinlon(ig)=sin(plon(ig))
1077	coslon(ig)=cos(plon(ig))
1078	ENDDO
1079	#endif
1080	! initialize variables in radinc_h
1081	call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
1082
1083	! allocate "comsoil_h" arrays
1084	call ini_comsoil_h(ngrid)
1085
1086	END SUBROUTINE inifis
1087
1088	END MODULE inifis_mod

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

Download in other formats:

Original Format