Context Navigation

inifis_mod.F90 @ 2283

Last change on this file since 2283 was 2283, checked in by jvatant, 5 years ago

Set the temperature boundaries and step for Planck function integration as input in callphys.def, for more flexibility.

+ User can now set them by tplanckmin, tplanckmax and dtplanck
+ Default values are a wide range 30-1500K
+ Add a sanity check in callcorrk instead of leaving out-of-bounds planckir indexes.

--JVO

File size: 27.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 mod_phys_lmdz_para, only : is_parallel
25
26	!=======================================================================
27	!
28	! purpose:
29	! -------
30	!
31	! Initialisation for the physical parametrisations of the LMD
32	! Generic Model.
33	!
34	! author: Frederic Hourdin 15 / 10 /93
35	! -------
36	! modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
37	! Ehouarn Millour (oct. 2008) tracers are now identified
38	! by their names and may not be contiguously
39	! stored in the q(:,:,:,:) array
40	! E.M. (june 2009) use getin routine to load parameters
41	!
42	!
43	! arguments:
44	! ----------
45	!
46	! input:
47	! ------
48	!
49	! ngrid Size of the horizontal grid.
50	! All internal loops are performed on that grid.
51	! nlayer Number of vertical layers.
52	! pdayref Day of reference for the simulation
53	! pday Number of days counted from the North. Spring
54	! equinoxe.
55	!
56	!=======================================================================
57	!
58	!-----------------------------------------------------------------------
59	! declarations:
60	! -------------
61	use datafile_mod, only: datadir
62	use ioipsl_getin_p_mod, only: getin_p
63	IMPLICIT NONE
64
65
66
67	REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep
68	INTEGER,INTENT(IN) :: nday
69	INTEGER,INTENT(IN) :: ngrid,nlayer,nq
70	REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid)
71	integer,intent(in) :: day_ini
72	INTEGER ig,ierr
73
74	EXTERNAL iniorbit,orbite
75	EXTERNAL SSUM
76	REAL SSUM
77
78	! Initialize flags lunout, prt_level, debug (in print_control_mod)
79	CALL init_print_control
80
81	! initialize constants in comcstfi_mod
82	rad=prad
83	cpp=pcpp
84	g=pg
85	r=pr
86	rcp=r/cpp
87	mugaz=8.314*1000./pr ! dummy init
88	pi=2.*asin(1.)
89	avocado = 6.02214179e23 ! added by RW
90
91	! Initialize some "temporal and calendar" related variables
92	#ifndef MESOSCALE
93	CALL init_time(day_ini,pdaysec,nday,ptimestep)
94	#endif
95
96	! read in some parameters from "run.def" for physics,
97	! or shared between dynamics and physics.
98	call getin_p("ecritphy",ecritphy) ! frequency of outputs in physics,
99	! in dynamical steps
100	call getin_p("day_step",day_step) ! number of dynamical steps per day
101	call getin_p("iphysiq",iphysiq) ! call physics every iphysiq dyn step
102
103	! do we read a startphy.nc file? (default: .true.)
104	call getin_p("startphy_file",startphy_file)
105
106	! --------------------------------------------------------------
107	! Reading the "callphys.def" file controlling some key options
108	! --------------------------------------------------------------
109
110	! check that 'callphys.def' file is around
111	OPEN(99,file='callphys.def',status='old',form='formatted',iostat=ierr)
112	CLOSE(99)
113	IF(ierr.EQ.0) iscallphys=.true. !iscallphys initialised as false in callkeys_mod module
114
115	!!! IF(ierr.EQ.0) THEN
116	IF(iscallphys) THEN
117	PRINT*
118	PRINT*
119	PRINT*,'--------------------------------------------'
120	PRINT*,' inifis: Parametres pour la physique (callphys.def)'
121	PRINT*,'--------------------------------------------'
122
123	write(,) "Directory where external input files are:"
124	! default 'datadir' is set in "datadir_mod"
125	call getin_p("datadir",datadir) ! default path
126	write(,) " datadir = ",trim(datadir)
127
128	write(,) "Run with or without tracer transport ?"
129	tracer=.false. ! default value
130	call getin_p("tracer",tracer)
131	write(,) " tracer = ",tracer
132
133	write(,) "Run with or without atm mass update ", &
134	" due to tracer evaporation/condensation?"
135	mass_redistrib=.false. ! default value
136	call getin_p("mass_redistrib",mass_redistrib)
137	write(,) " mass_redistrib = ",mass_redistrib
138
139	write(,) "Diurnal cycle ?"
140	write(,) "(if diurnal=false, diurnal averaged solar heating)"
141	diurnal=.true. ! default value
142	call getin_p("diurnal",diurnal)
143	write(,) " diurnal = ",diurnal
144
145	write(,) "Seasonal cycle ?"
146	write(,) "(if season=false, Ls stays constant, to value ", &
147	"set in 'start'"
148	season=.true. ! default value
149	call getin_p("season",season)
150	write(,) " season = ",season
151
152	write(,) "Tidally resonant rotation ?"
153	tlocked=.false. ! default value
154	call getin_p("tlocked",tlocked)
155	write(,) "tlocked = ",tlocked
156
157	write(,) "Saturn ring shadowing ?"
158	rings_shadow = .false.
159	call getin_p("rings_shadow", rings_shadow)
160	write(,) "rings_shadow = ", rings_shadow
161
162	write(,) "Compute latitude-dependent gravity field?"
163	oblate = .false.
164	call getin_p("oblate", oblate)
165	write(,) "oblate = ", oblate
166
167	write(,) "Flattening of the planet (a-b)/a "
168	flatten = 0.0
169	call getin_p("flatten", flatten)
170	write(,) "flatten = ", flatten
171
172
173	write(,) "Needed if oblate=.true.: J2"
174	J2 = 0.0
175	call getin_p("J2", J2)
176	write(,) "J2 = ", J2
177
178	write(,) "Needed if oblate=.true.: Planet mass (*1e24 kg)"
179	MassPlanet = 0.0
180	call getin_p("MassPlanet", MassPlanet)
181	write(,) "MassPlanet = ", MassPlanet
182
183	write(,) "Needed if oblate=.true.: Planet mean radius (m)"
184	Rmean = 0.0
185	call getin_p("Rmean", Rmean)
186	write(,) "Rmean = ", Rmean
187
188	! Test of incompatibility:
189	! if tlocked, then diurnal should be false
190	if (tlocked.and.diurnal) then
191	print*,'If diurnal=true, we should turn off tlocked.'
192	stop
193	endif
194
195	write(,) "Tidal resonance ratio ?"
196	nres=0 ! default value
197	call getin_p("nres",nres)
198	write(,) "nres = ",nres
199
200	write(,) "Write some extra output to the screen ?"
201	lwrite=.false. ! default value
202	call getin_p("lwrite",lwrite)
203	write(,) " lwrite = ",lwrite
204
205	write(,) "Save statistics in file stats.nc ?"
206	callstats=.true. ! default value
207	call getin_p("callstats",callstats)
208	write(,) " callstats = ",callstats
209
210	write(,) "Test energy conservation of model physics ?"
211	enertest=.false. ! default value
212	call getin_p("enertest",enertest)
213	write(,) " enertest = ",enertest
214
215	write(,) "Check to see if cpp values used match gases.def ?"
216	check_cpp_match=.true. ! default value
217	call getin_p("check_cpp_match",check_cpp_match)
218	write(,) " check_cpp_match = ",check_cpp_match
219
220	write(,) "call radiative transfer ?"
221	callrad=.true. ! default value
222	call getin_p("callrad",callrad)
223	write(,) " callrad = ",callrad
224
225	write(,) "call correlated-k radiative transfer ?"
226	corrk=.true. ! default value
227	call getin_p("corrk",corrk)
228	write(,) " corrk = ",corrk
229
230	write(,) "prohibit calculations outside corrk T grid?"
231	strictboundcorrk=.true. ! default value
232	call getin_p("strictboundcorrk",strictboundcorrk)
233	write(,) "strictboundcorrk = ",strictboundcorrk
234
235	write(,) "Minimum atmospheric temperature for Planck function integration ?"
236	tplanckmin=30.0 ! default value
237	call getin_p("tplanckmin",tplanckmin)
238	write(,) " tplanckmin = ",tplanckmin
239
240	write(,) "Maximum atmospheric temperature for Planck function integration ?"
241	tplanckmax=1500.0 ! default value
242	call getin_p("tplanckmax",tplanckmax)
243	write(,) " tplanckmax = ",tplanckmax
244
245	write(,) "Temperature step for Planck function integration ?"
246	dtplanck=0.1 ! default value
247	call getin_p("dtplanck",dtplanck)
248	write(,) " dtplanck = ",dtplanck
249
250	write(,) "call gaseous absorption in the visible bands?", &
251	"(matters only if callrad=T)"
252	callgasvis=.false. ! default value
253	call getin_p("callgasvis",callgasvis)
254	write(,) " callgasvis = ",callgasvis
255
256	write(,) "call continuum opacities in radiative transfer ?", &
257	"(matters only if callrad=T)"
258	continuum=.true. ! default value
259	call getin_p("continuum",continuum)
260	write(,) " continuum = ",continuum
261
262	write(,) "use analytic function for H2O continuum ?"
263	H2Ocont_simple=.false. ! default value
264	call getin_p("H2Ocont_simple",H2Ocont_simple)
265	write(,) " H2Ocont_simple = ",H2Ocont_simple
266
267	write(,) "version for H2H2 CIA file ?"
268	versH2H2cia=2011 ! default value (should be 2018 but retrocompatibility first)
269	call getin_p("versH2H2cia",versH2H2cia)
270	write(,) " versH2H2cia = ",versH2H2cia
271	! Sanity check
272	if (versH2H2cia.ne.2011 .and. versH2H2cia.ne.2018) then
273	print*,'Error: Choose a correct value (2011 or 2018) for versH2H2cia !'
274	call abort
275	endif
276
277	write(,) "call turbulent vertical diffusion ?"
278	calldifv=.true. ! default value
279	call getin_p("calldifv",calldifv)
280	write(,) " calldifv = ",calldifv
281
282	write(,) "use turbdiff instead of vdifc ?"
283	UseTurbDiff=.true. ! default value
284	call getin_p("UseTurbDiff",UseTurbDiff)
285	write(,) " UseTurbDiff = ",UseTurbDiff
286
287	write(,) "call convective adjustment ?"
288	calladj=.true. ! default value
289	call getin_p("calladj",calladj)
290	write(,) " calladj = ",calladj
291
292	write(,) "call thermal plume model ?"
293	calltherm=.false. ! default value
294	call getin_p("calltherm",calltherm)
295	write(,) " calltherm = ",calltherm
296
297	write(,) "call CO2 condensation ?"
298	co2cond=.false. ! default value
299	call getin_p("co2cond",co2cond)
300	write(,) " co2cond = ",co2cond
301	! Test of incompatibility
302	if (co2cond.and.(.not.tracer)) then
303	print*,'We need a CO2 ice tracer to condense CO2'
304	call abort
305	endif
306
307	write(,) "CO2 supersaturation level ?"
308	co2supsat=1.0 ! default value
309	call getin_p("co2supsat",co2supsat)
310	write(,) " co2supsat = ",co2supsat
311
312	write(,) "Radiative timescale for Newtonian cooling ?"
313	tau_relax=30. ! default value
314	call getin_p("tau_relax",tau_relax)
315	write(,) " tau_relax = ",tau_relax
316	tau_relax=tau_relax243600 ! convert Earth days --> seconds
317
318	write(,)"call thermal conduction in the soil ?"
319	callsoil=.true. ! default value
320	call getin_p("callsoil",callsoil)
321	write(,) " callsoil = ",callsoil
322
323	write(,)"Rad transfer is computed every iradia", &
324	" physical timestep"
325	iradia=1 ! default value
326	call getin_p("iradia",iradia)
327	write(,)" iradia = ",iradia
328
329	write(,)"Rayleigh scattering ?"
330	rayleigh=.false.
331	call getin_p("rayleigh",rayleigh)
332	write(,)" rayleigh = ",rayleigh
333
334	write(,) "Use blackbody for stellar spectrum ?"
335	stelbbody=.false. ! default value
336	call getin_p("stelbbody",stelbbody)
337	write(,) " stelbbody = ",stelbbody
338
339	write(,) "Stellar blackbody temperature ?"
340	stelTbb=5800.0 ! default value
341	call getin_p("stelTbb",stelTbb)
342	write(,) " stelTbb = ",stelTbb
343
344	write(,)"Output mean OLR in 1D?"
345	meanOLR=.false.
346	call getin_p("meanOLR",meanOLR)
347	write(,)" meanOLR = ",meanOLR
348
349	write(,)"Output spectral OLR in 3D?"
350	specOLR=.false.
351	call getin_p("specOLR",specOLR)
352	write(,)" specOLR = ",specOLR
353
354	write(,)"Output diagnostic optical thickness attenuation exp(-tau)?"
355	diagdtau=.false.
356	call getin_p("diagdtau",diagdtau)
357	write(,)" diagdtau = ",diagdtau
358
359	write(,)"Operate in kastprof mode?"
360	kastprof=.false.
361	call getin_p("kastprof",kastprof)
362	write(,)" kastprof = ",kastprof
363
364	write(,)"Uniform absorption in radiative transfer?"
365	graybody=.false.
366	call getin_p("graybody",graybody)
367	write(,)" graybody = ",graybody
368
369	! Soil model
370	write(,)"Number of sub-surface layers for soil scheme?"
371	! default value of nsoilmx set in comsoil_h
372	call getin_p("nsoilmx",nsoilmx)
373	write(,)" nsoilmx=",nsoilmx
374
375	write(,)"Thickness of topmost soil layer (m)?"
376	! default value of lay1_soil set in comsoil_h
377	call getin_p("lay1_soil",lay1_soil)
378	write(,)" lay1_soil = ",lay1_soil
379
380	write(,)"Coefficient for soil layer thickness distribution?"
381	! default value of alpha_soil set in comsoil_h
382	call getin_p("alpha_soil",alpha_soil)
383	write(,)" alpha_soil = ",alpha_soil
384
385	! Slab Ocean
386	write(,) "Use slab-ocean ?"
387	ok_slab_ocean=.false. ! default value
388	call getin_p("ok_slab_ocean",ok_slab_ocean)
389	write(,) "ok_slab_ocean = ",ok_slab_ocean
390	! Sanity check: for now slab oncean only works in serial mode
391	if (ok_slab_ocean.and.is_parallel) then
392	write(,) " Error: slab ocean should only be used in serial mode!"
393	call abort
394	endif
395
396	write(,) "Use slab-sea-ice ?"
397	ok_slab_sic=.true. ! default value
398	call getin_p("ok_slab_sic",ok_slab_sic)
399	write(,) "ok_slab_sic = ",ok_slab_sic
400
401	write(,) "Use heat transport for the ocean ?"
402	ok_slab_heat_transp=.true. ! default value
403	call getin_p("ok_slab_heat_transp",ok_slab_heat_transp)
404	write(,) "ok_slab_heat_transp = ",ok_slab_heat_transp
405
406	! Photochemistry and chemistry in the thermosphere
407
408	write(,) "Use photochemistry ?"
409	photochem=.false. ! default value
410	call getin_p("photochem",photochem)
411	write(,) "photochem = ",photochem
412
413	write(,)"Production of haze ?"
414	haze=.false. ! default value
415	call getin_p("haze",haze)
416	write(,)" haze = ",haze
417
418
419	! Test of incompatibility:
420	! if kastprof used, we must be in 1D
421	if (kastprof.and.(ngrid.gt.1)) then
422	print*,'kastprof can only be used in 1D!'
423	call abort
424	endif
425
426	write(,)"Stratospheric temperature for kastprof mode?"
427	Tstrat=167.0
428	call getin_p("Tstrat",Tstrat)
429	write(,)" Tstrat = ",Tstrat
430
431	write(,)"Remove lower boundary?"
432	nosurf=.false.
433	call getin_p("nosurf",nosurf)
434	write(,)" nosurf = ",nosurf
435
436	! Tests of incompatibility:
437	if (nosurf.and.callsoil) then
438	print*,'nosurf not compatible with soil scheme!'
439	print*,'... got to make a choice!'
440	call abort
441	endif
442
443	write(,)"Add an internal heat flux?", &
444	"... matters only if callsoil=F"
445	intheat=0.
446	call getin_p("intheat",intheat)
447	write(,)" intheat = ",intheat
448
449	write(,)"Use Newtonian cooling for radiative transfer?"
450	newtonian=.false.
451	call getin_p("newtonian",newtonian)
452	write(,)" newtonian = ",newtonian
453
454	! Tests of incompatibility:
455	if (newtonian.and.corrk) then
456	print*,'newtonian not compatible with correlated-k!'
457	call abort
458	endif
459	if (newtonian.and.calladj) then
460	print*,'newtonian not compatible with adjustment!'
461	call abort
462	endif
463	if (newtonian.and.calldifv) then
464	print*,'newtonian not compatible with a boundary layer!'
465	call abort
466	endif
467
468	write(,)"Test physics timescale in 1D?"
469	testradtimes=.false.
470	call getin_p("testradtimes",testradtimes)
471	write(,)" testradtimes = ",testradtimes
472
473	! Test of incompatibility:
474	! if testradtimes used, we must be in 1D
475	if (testradtimes.and.(ngrid.gt.1)) then
476	print*,'testradtimes can only be used in 1D!'
477	call abort
478	endif
479
480	write(,)"Default planetary temperature?"
481	tplanet=215.0
482	call getin_p("tplanet",tplanet)
483	write(,)" tplanet = ",tplanet
484
485	write(,)"Which star?"
486	startype=1 ! default value = Sol
487	call getin_p("startype",startype)
488	write(,)" startype = ",startype
489
490	write(,)"Value of stellar flux at 1 AU?"
491	Fat1AU=1356.0 ! default value = Sol today
492	call getin_p("Fat1AU",Fat1AU)
493	write(,)" Fat1AU = ",Fat1AU
494
495
496	! TRACERS:
497
498	write(,)"Varying H2O cloud fraction?"
499	CLFvarying=.false. ! default value
500	call getin_p("CLFvarying",CLFvarying)
501	write(,)" CLFvarying = ",CLFvarying
502
503	write(,)"Value of fixed H2O cloud fraction?"
504	CLFfixval=1.0 ! default value
505	call getin_p("CLFfixval",CLFfixval)
506	write(,)" CLFfixval = ",CLFfixval
507
508	write(,)"fixed radii for Cloud particles?"
509	radfixed=.false. ! default value
510	call getin_p("radfixed",radfixed)
511	write(,)" radfixed = ",radfixed
512
513	if(kastprof)then
514	radfixed=.true.
515	endif
516
517	write(,)"Number mixing ratio of CO2 ice particles:"
518	Nmix_co2=1.e6 ! default value
519	call getin_p("Nmix_co2",Nmix_co2)
520	write(,)" Nmix_co2 = ",Nmix_co2
521
522	! write(,)"Number of radiatively active aerosols:"
523	! naerkind=0. ! default value
524	! call getin_p("naerkind",naerkind)
525	! write(,)" naerkind = ",naerkind
526
527	write(,)"Opacity of dust (if used):"
528	dusttau=0. ! default value
529	call getin_p("dusttau",dusttau)
530	write(,)" dusttau = ",dusttau
531
532	write(,)"Radiatively active CO2 aerosols?"
533	aeroco2=.false. ! default value
534	call getin_p("aeroco2",aeroco2)
535	write(,)" aeroco2 = ",aeroco2
536
537	write(,)"Fixed CO2 aerosol distribution?"
538	aerofixco2=.false. ! default value
539	call getin_p("aerofixco2",aerofixco2)
540	write(,)" aerofixco2 = ",aerofixco2
541
542	write(,)"Radiatively active water ice?"
543	aeroh2o=.false. ! default value
544	call getin_p("aeroh2o",aeroh2o)
545	write(,)" aeroh2o = ",aeroh2o
546
547	write(,)"Fixed H2O aerosol distribution?"
548	aerofixh2o=.false. ! default value
549	call getin_p("aerofixh2o",aerofixh2o)
550	write(,)" aerofixh2o = ",aerofixh2o
551
552	write(,)"Radiatively active sulfuric acid aerosols?"
553	aeroh2so4=.false. ! default value
554	call getin_p("aeroh2so4",aeroh2so4)
555	write(,)" aeroh2so4 = ",aeroh2so4
556
557	!=================================
558
559	write(,)"Radiatively active two-layer aerosols?"
560	aeroback2lay=.false. ! default value
561	call getin_p("aeroback2lay",aeroback2lay)
562	write(,)" aeroback2lay = ",aeroback2lay
563
564	write(,)"Radiatively active ammonia cloud?"
565	aeronh3=.false. ! default value
566	call getin_p("aeronh3",aeronh3)
567	write(,)" aeronh3 = ",aeronh3
568
569	write(,)"Radiatively active auroral aerosols?"
570	aeroaurora=.false. ! default value
571	call getin_p("aeroaurora",aeroaurora)
572	write(,)" aeroaurora = ",aeroaurora
573
574	write(,)"TWOLAY AEROSOL: total optical depth ", &
575	"in the tropospheric layer (visible)"
576	obs_tau_col_tropo=8.D0
577	call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
578	write(,)" obs_tau_col_tropo = ",obs_tau_col_tropo
579
580	write(,)"TWOLAY AEROSOL: total optical depth ", &
581	"in the stratospheric layer (visible)"
582	obs_tau_col_strato=0.08D0
583	call getin_p("obs_tau_col_strato",obs_tau_col_strato)
584	write(,)" obs_tau_col_strato = ",obs_tau_col_strato
585
586	write(,)"TWOLAY AEROSOL: optprop_back2lay_vis?"
587	optprop_back2lay_vis = 'optprop_saturn_vis_n20.dat'
588	call getin_p("optprop_back2lay_vis",optprop_back2lay_vis)
589	write(,)" optprop_back2lay_vis = ",optprop_back2lay_vis
590
591	write(,)"TWOLAY AEROSOL: optprop_back2lay_ir?"
592	optprop_back2lay_ir = 'optprop_saturn_ir_n20.dat'
593	call getin_p("optprop_back2lay_ir",optprop_back2lay_ir)
594	write(,)" optprop_back2lay_ir = ",optprop_back2lay_ir
595
596	write(,)"TWOLAY AEROSOL: pres_bottom_tropo? in pa"
597	pres_bottom_tropo=66000.0
598	call getin_p("pres_bottom_tropo",pres_bottom_tropo)
599	write(,)" pres_bottom_tropo = ",pres_bottom_tropo
600
601	write(,)"TWOLAY AEROSOL: pres_top_tropo? in pa"
602	pres_top_tropo=18000.0
603	call getin_p("pres_top_tropo",pres_top_tropo)
604	write(,)" pres_top_tropo = ",pres_top_tropo
605
606	write(,)"TWOLAY AEROSOL: pres_bottom_strato? in pa"
607	pres_bottom_strato=2000.0
608	call getin_p("pres_bottom_strato",pres_bottom_strato)
609	write(,)" pres_bottom_strato = ",pres_bottom_strato
610
611	! Sanity check
612	if (pres_bottom_strato .gt. pres_top_tropo) then
613	print*,'Error : TWOLAY AEROSOL, Please ensure that in callphys.def'
614	print*,'you have pres_top_tropo > pres_bottom_strato !'
615	stop
616	endif
617
618	write(,)"TWOLAY AEROSOL: pres_top_strato? in pa"
619	pres_top_strato=100.0
620	call getin_p("pres_top_strato",pres_top_strato)
621	write(,)" pres_top_strato = ",pres_top_strato
622
623	write(,)"TWOLAY AEROSOL: particle size in the ", &
624	"tropospheric layer, in meters"
625	size_tropo=2.e-6
626	call getin_p("size_tropo",size_tropo)
627	write(,)" size_tropo = ",size_tropo
628
629	write(,)"TWOLAY AEROSOL: particle size in the ", &
630	"stratospheric layer, in meters"
631	size_strato=1.e-7
632	call getin_p("size_strato",size_strato)
633	write(,)" size_strato = ",size_strato
634
635	write(,)"NH3 (thin) cloud: total optical depth"
636	tau_nh3_cloud=7.D0
637	call getin_p("tau_nh3_cloud",tau_nh3_cloud)
638	write(,)" tau_nh3_cloud = ",tau_nh3_cloud
639
640	write(,)"NH3 (thin) cloud pressure level"
641	pres_nh3_cloud=7.D0
642	call getin_p("pres_nh3_cloud",pres_nh3_cloud)
643	write(,)" pres_nh3_cloud = ",pres_nh3_cloud
644
645	write(,)"NH3 (thin) cloud: particle sizes"
646	size_nh3_cloud=3.e-6
647	call getin_p("size_nh3_cloud",size_nh3_cloud)
648	write(,)" size_nh3_cloud = ",size_nh3_cloud
649
650	!=================================
651
652	write(,)"Cloud pressure level (with kastprof only):"
653	cloudlvl=0. ! default value
654	call getin_p("cloudlvl",cloudlvl)
655	write(,)" cloudlvl = ",cloudlvl
656
657	write(,)"Is the variable gas species radiatively active?"
658	Tstrat=167.0
659	varactive=.false.
660	call getin_p("varactive",varactive)
661	write(,)" varactive = ",varactive
662
663	write(,)"Is the variable gas species distribution set?"
664	varfixed=.false.
665	call getin_p("varfixed",varfixed)
666	write(,)" varfixed = ",varfixed
667
668	write(,)"What is the saturation % of the variable species?"
669	satval=0.8
670	call getin_p("satval",satval)
671	write(,)" satval = ",satval
672
673
674	! Test of incompatibility:
675	! if varactive, then varfixed should be false
676	if (varactive.and.varfixed) then
677	print*,'if varactive, varfixed must be OFF!'
678	stop
679	endif
680
681	write(,) "Gravitationnal sedimentation ?"
682	sedimentation=.false. ! default value
683	call getin_p("sedimentation",sedimentation)
684	write(,) " sedimentation = ",sedimentation
685
686	write(,) "Compute water cycle ?"
687	water=.false. ! default value
688	call getin_p("water",water)
689	write(,) " water = ",water
690
691	! Test of incompatibility:
692	! if water is true, there should be at least a tracer
693	if (water.and.(.not.tracer)) then
694	print*,'if water is ON, tracer must be ON too!'
695	stop
696	endif
697
698	write(,) "Include water condensation ?"
699	watercond=.false. ! default value
700	call getin_p("watercond",watercond)
701	write(,) " watercond = ",watercond
702
703	! Test of incompatibility:
704	! if watercond is used, then water should be used too
705	if (watercond.and.(.not.water)) then
706	print*,'if watercond is used, water should be used too'
707	stop
708	endif
709
710	write(,) "Include water precipitation ?"
711	waterrain=.false. ! default value
712	call getin_p("waterrain",waterrain)
713	write(,) " waterrain = ",waterrain
714
715	write(,) "Include surface hydrology ?"
716	hydrology=.false. ! default value
717	call getin_p("hydrology",hydrology)
718	write(,) " hydrology = ",hydrology
719
720	write(,) "Evolve surface water sources ?"
721	sourceevol=.false. ! default value
722	call getin_p("sourceevol",sourceevol)
723	write(,) " sourceevol = ",sourceevol
724
725	write(,) "Ice evolution timestep ?"
726	icetstep=100.0 ! default value
727	call getin_p("icetstep",icetstep)
728	write(,) " icetstep = ",icetstep
729
730	write(,) "Spectral Dependant albedo ?"
731	albedo_spectral_mode=.false. ! default value
732	call getin_p("albedo_spectral_mode",albedo_spectral_mode)
733	write(,) " albedo_spectral_mode = ",albedo_spectral_mode
734
735	write(,) "Snow albedo ?"
736	write(,) "If albedo_spectral_mode=.true., then this "
737	write(,) "corresponds to the 0.5 microns snow albedo."
738	albedosnow=0.5 ! default value
739	call getin_p("albedosnow",albedosnow)
740	write(,) " albedosnow = ",albedosnow
741
742	write(,) "CO2 ice albedo ?"
743	albedoco2ice=0.5 ! default value
744	call getin_p("albedoco2ice",albedoco2ice)
745	write(,) " albedoco2ice = ",albedoco2ice
746
747	write(,) "Maximum ice thickness ?"
748	maxicethick=2.0 ! default value
749	call getin_p("maxicethick",maxicethick)
750	write(,) " maxicethick = ",maxicethick
751
752	write(,) "Freezing point of seawater ?"
753	Tsaldiff=-1.8 ! default value
754	call getin_p("Tsaldiff",Tsaldiff)
755	write(,) " Tsaldiff = ",Tsaldiff
756
757	write(,) "Does user want to force cpp and mugaz?"
758	force_cpp=.false. ! default value
759	call getin_p("force_cpp",force_cpp)
760	write(,) " force_cpp = ",force_cpp
761
762	if (force_cpp) then
763	mugaz = -99999.
764	PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
765	call getin_p("mugaz",mugaz)
766	IF (mugaz.eq.-99999.) THEN
767	PRINT *, "mugaz must be set if force_cpp = T"
768	STOP
769	ELSE
770	write(,) "mugaz=",mugaz
771	ENDIF
772	cpp = -99999.
773	PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
774	call getin_p("cpp",cpp)
775	IF (cpp.eq.-99999.) THEN
776	PRINT *, "cpp must be set if force_cpp = T"
777	STOP
778	ELSE
779	write(,) "cpp=",cpp
780	ENDIF
781	endif ! of if (force_cpp)
782	call su_gases
783	call calc_cpp_mugaz
784
785	PRINT*,'--------------------------------------------'
786	PRINT*
787	PRINT*
788	ELSE
789	write(,)
790	write(,) 'Cannot read file callphys.def. Is it here ?'
791	stop
792	ENDIF ! of IF(iscallphys)
793
794	PRINT*
795	PRINT*,'inifis: daysec',daysec
796	PRINT*
797	PRINT*,'inifis: The radiative transfer is computed:'
798	PRINT*,' each ',iradia,' physical time-step'
799	PRINT,' or each ',iradiadtphys,' seconds'
800	PRINT*
801
802
803	!-----------------------------------------------------------------------
804	! Some more initialization:
805	! ------------------------
806
807	! Initializations for comgeomfi_h
808	#ifndef MESOSCALE
809	totarea=SSUM(ngrid,parea,1)
810	call planetwide_sumval(parea,totarea_planet)
811
812	!! those are defined in comdiurn_h.F90
813	IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
814	IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
815	IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
816	IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
817
818	DO ig=1,ngrid
819	sinlat(ig)=sin(plat(ig))
820	coslat(ig)=cos(plat(ig))
821	sinlon(ig)=sin(plon(ig))
822	coslon(ig)=cos(plon(ig))
823	ENDDO
824	#endif
825	! initialize variables in radinc_h
826	call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
827
828	! allocate "comsoil_h" arrays
829	call ini_comsoil_h(ngrid)
830
831	END SUBROUTINE inifis
832
833	END MODULE inifis_mod

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

Download in other formats:

Original Format