Context Navigation

inifis_mod.F90 @ 2254

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

Fix a bug in the 2-layers aerosol scheme in aeropacity.F90
-> It appeared that the normalization for total column optical depth was a

wrong way to go, the actual strato and tropo values weren't matcihing the
input values, and one of the main consequence was to strongly underestimate
the actual strato optical depth compared to input value in most cases.

-> Now we ensure the normalization to the input values for both layers, while

keeping a dependence in dP for each GCM layer.

-> Also add a sanity check ensuring pres_top_tropo > pres_bottom_strato.
--JVO

File size: 26.6 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(,) "call gaseous absorption in the visible bands?", &
236	"(matters only if callrad=T)"
237	callgasvis=.false. ! default value
238	call getin_p("callgasvis",callgasvis)
239	write(,) " callgasvis = ",callgasvis
240
241	write(,) "call continuum opacities in radiative transfer ?", &
242	"(matters only if callrad=T)"
243	continuum=.true. ! default value
244	call getin_p("continuum",continuum)
245	write(,) " continuum = ",continuum
246
247	write(,) "use analytic function for H2O continuum ?"
248	H2Ocont_simple=.false. ! default value
249	call getin_p("H2Ocont_simple",H2Ocont_simple)
250	write(,) " H2Ocont_simple = ",H2Ocont_simple
251
252	write(,) "version for H2H2 CIA file ?"
253	versH2H2cia=2011 ! default value (should be 2018 but retrocompatibility first)
254	call getin_p("versH2H2cia",versH2H2cia)
255	write(,) " versH2H2cia = ",versH2H2cia
256	! Sanity check
257	if (versH2H2cia.ne.2011 .and. versH2H2cia.ne.2018) then
258	print*,'Error: Choose a correct value (2011 or 2018) for versH2H2cia !'
259	call abort
260	endif
261
262	write(,) "call turbulent vertical diffusion ?"
263	calldifv=.true. ! default value
264	call getin_p("calldifv",calldifv)
265	write(,) " calldifv = ",calldifv
266
267	write(,) "use turbdiff instead of vdifc ?"
268	UseTurbDiff=.true. ! default value
269	call getin_p("UseTurbDiff",UseTurbDiff)
270	write(,) " UseTurbDiff = ",UseTurbDiff
271
272	write(,) "call convective adjustment ?"
273	calladj=.true. ! default value
274	call getin_p("calladj",calladj)
275	write(,) " calladj = ",calladj
276
277	write(,) "call thermal plume model ?"
278	calltherm=.false. ! default value
279	call getin_p("calltherm",calltherm)
280	write(,) " calltherm = ",calltherm
281
282	write(,) "call CO2 condensation ?"
283	co2cond=.false. ! default value
284	call getin_p("co2cond",co2cond)
285	write(,) " co2cond = ",co2cond
286	! Test of incompatibility
287	if (co2cond.and.(.not.tracer)) then
288	print*,'We need a CO2 ice tracer to condense CO2'
289	call abort
290	endif
291
292	write(,) "CO2 supersaturation level ?"
293	co2supsat=1.0 ! default value
294	call getin_p("co2supsat",co2supsat)
295	write(,) " co2supsat = ",co2supsat
296
297	write(,) "Radiative timescale for Newtonian cooling ?"
298	tau_relax=30. ! default value
299	call getin_p("tau_relax",tau_relax)
300	write(,) " tau_relax = ",tau_relax
301	tau_relax=tau_relax243600 ! convert Earth days --> seconds
302
303	write(,)"call thermal conduction in the soil ?"
304	callsoil=.true. ! default value
305	call getin_p("callsoil",callsoil)
306	write(,) " callsoil = ",callsoil
307
308	write(,)"Rad transfer is computed every iradia", &
309	" physical timestep"
310	iradia=1 ! default value
311	call getin_p("iradia",iradia)
312	write(,)" iradia = ",iradia
313
314	write(,)"Rayleigh scattering ?"
315	rayleigh=.false.
316	call getin_p("rayleigh",rayleigh)
317	write(,)" rayleigh = ",rayleigh
318
319	write(,) "Use blackbody for stellar spectrum ?"
320	stelbbody=.false. ! default value
321	call getin_p("stelbbody",stelbbody)
322	write(,) " stelbbody = ",stelbbody
323
324	write(,) "Stellar blackbody temperature ?"
325	stelTbb=5800.0 ! default value
326	call getin_p("stelTbb",stelTbb)
327	write(,) " stelTbb = ",stelTbb
328
329	write(,)"Output mean OLR in 1D?"
330	meanOLR=.false.
331	call getin_p("meanOLR",meanOLR)
332	write(,)" meanOLR = ",meanOLR
333
334	write(,)"Output spectral OLR in 3D?"
335	specOLR=.false.
336	call getin_p("specOLR",specOLR)
337	write(,)" specOLR = ",specOLR
338
339	write(,)"Output diagnostic optical thickness attenuation exp(-tau)?"
340	diagdtau=.false.
341	call getin_p("diagdtau",diagdtau)
342	write(,)" diagdtau = ",diagdtau
343
344	write(,)"Operate in kastprof mode?"
345	kastprof=.false.
346	call getin_p("kastprof",kastprof)
347	write(,)" kastprof = ",kastprof
348
349	write(,)"Uniform absorption in radiative transfer?"
350	graybody=.false.
351	call getin_p("graybody",graybody)
352	write(,)" graybody = ",graybody
353
354	! Soil model
355	write(,)"Number of sub-surface layers for soil scheme?"
356	! default value of nsoilmx set in comsoil_h
357	call getin_p("nsoilmx",nsoilmx)
358	write(,)" nsoilmx=",nsoilmx
359
360	write(,)"Thickness of topmost soil layer (m)?"
361	! default value of lay1_soil set in comsoil_h
362	call getin_p("lay1_soil",lay1_soil)
363	write(,)" lay1_soil = ",lay1_soil
364
365	write(,)"Coefficient for soil layer thickness distribution?"
366	! default value of alpha_soil set in comsoil_h
367	call getin_p("alpha_soil",alpha_soil)
368	write(,)" alpha_soil = ",alpha_soil
369
370	! Slab Ocean
371	write(,) "Use slab-ocean ?"
372	ok_slab_ocean=.false. ! default value
373	call getin_p("ok_slab_ocean",ok_slab_ocean)
374	write(,) "ok_slab_ocean = ",ok_slab_ocean
375	! Sanity check: for now slab oncean only works in serial mode
376	if (ok_slab_ocean.and.is_parallel) then
377	write(,) " Error: slab ocean should only be used in serial mode!"
378	call abort
379	endif
380
381	write(,) "Use slab-sea-ice ?"
382	ok_slab_sic=.true. ! default value
383	call getin_p("ok_slab_sic",ok_slab_sic)
384	write(,) "ok_slab_sic = ",ok_slab_sic
385
386	write(,) "Use heat transport for the ocean ?"
387	ok_slab_heat_transp=.true. ! default value
388	call getin_p("ok_slab_heat_transp",ok_slab_heat_transp)
389	write(,) "ok_slab_heat_transp = ",ok_slab_heat_transp
390
391	! Photochemistry and chemistry in the thermosphere
392
393	write(,) "Use photochemistry ?"
394	photochem=.false. ! default value
395	call getin_p("photochem",photochem)
396	write(,) "photochem = ",photochem
397
398	write(,)"Production of haze ?"
399	haze=.false. ! default value
400	call getin_p("haze",haze)
401	write(,)" haze = ",haze
402
403
404	! Test of incompatibility:
405	! if kastprof used, we must be in 1D
406	if (kastprof.and.(ngrid.gt.1)) then
407	print*,'kastprof can only be used in 1D!'
408	call abort
409	endif
410
411	write(,)"Stratospheric temperature for kastprof mode?"
412	Tstrat=167.0
413	call getin_p("Tstrat",Tstrat)
414	write(,)" Tstrat = ",Tstrat
415
416	write(,)"Remove lower boundary?"
417	nosurf=.false.
418	call getin_p("nosurf",nosurf)
419	write(,)" nosurf = ",nosurf
420
421	! Tests of incompatibility:
422	if (nosurf.and.callsoil) then
423	print*,'nosurf not compatible with soil scheme!'
424	print*,'... got to make a choice!'
425	call abort
426	endif
427
428	write(,)"Add an internal heat flux?", &
429	"... matters only if callsoil=F"
430	intheat=0.
431	call getin_p("intheat",intheat)
432	write(,)" intheat = ",intheat
433
434	write(,)"Use Newtonian cooling for radiative transfer?"
435	newtonian=.false.
436	call getin_p("newtonian",newtonian)
437	write(,)" newtonian = ",newtonian
438
439	! Tests of incompatibility:
440	if (newtonian.and.corrk) then
441	print*,'newtonian not compatible with correlated-k!'
442	call abort
443	endif
444	if (newtonian.and.calladj) then
445	print*,'newtonian not compatible with adjustment!'
446	call abort
447	endif
448	if (newtonian.and.calldifv) then
449	print*,'newtonian not compatible with a boundary layer!'
450	call abort
451	endif
452
453	write(,)"Test physics timescale in 1D?"
454	testradtimes=.false.
455	call getin_p("testradtimes",testradtimes)
456	write(,)" testradtimes = ",testradtimes
457
458	! Test of incompatibility:
459	! if testradtimes used, we must be in 1D
460	if (testradtimes.and.(ngrid.gt.1)) then
461	print*,'testradtimes can only be used in 1D!'
462	call abort
463	endif
464
465	write(,)"Default planetary temperature?"
466	tplanet=215.0
467	call getin_p("tplanet",tplanet)
468	write(,)" tplanet = ",tplanet
469
470	write(,)"Which star?"
471	startype=1 ! default value = Sol
472	call getin_p("startype",startype)
473	write(,)" startype = ",startype
474
475	write(,)"Value of stellar flux at 1 AU?"
476	Fat1AU=1356.0 ! default value = Sol today
477	call getin_p("Fat1AU",Fat1AU)
478	write(,)" Fat1AU = ",Fat1AU
479
480
481	! TRACERS:
482
483	write(,)"Varying H2O cloud fraction?"
484	CLFvarying=.false. ! default value
485	call getin_p("CLFvarying",CLFvarying)
486	write(,)" CLFvarying = ",CLFvarying
487
488	write(,)"Value of fixed H2O cloud fraction?"
489	CLFfixval=1.0 ! default value
490	call getin_p("CLFfixval",CLFfixval)
491	write(,)" CLFfixval = ",CLFfixval
492
493	write(,)"fixed radii for Cloud particles?"
494	radfixed=.false. ! default value
495	call getin_p("radfixed",radfixed)
496	write(,)" radfixed = ",radfixed
497
498	if(kastprof)then
499	radfixed=.true.
500	endif
501
502	write(,)"Number mixing ratio of CO2 ice particles:"
503	Nmix_co2=1.e6 ! default value
504	call getin_p("Nmix_co2",Nmix_co2)
505	write(,)" Nmix_co2 = ",Nmix_co2
506
507	! write(,)"Number of radiatively active aerosols:"
508	! naerkind=0. ! default value
509	! call getin_p("naerkind",naerkind)
510	! write(,)" naerkind = ",naerkind
511
512	write(,)"Opacity of dust (if used):"
513	dusttau=0. ! default value
514	call getin_p("dusttau",dusttau)
515	write(,)" dusttau = ",dusttau
516
517	write(,)"Radiatively active CO2 aerosols?"
518	aeroco2=.false. ! default value
519	call getin_p("aeroco2",aeroco2)
520	write(,)" aeroco2 = ",aeroco2
521
522	write(,)"Fixed CO2 aerosol distribution?"
523	aerofixco2=.false. ! default value
524	call getin_p("aerofixco2",aerofixco2)
525	write(,)" aerofixco2 = ",aerofixco2
526
527	write(,)"Radiatively active water ice?"
528	aeroh2o=.false. ! default value
529	call getin_p("aeroh2o",aeroh2o)
530	write(,)" aeroh2o = ",aeroh2o
531
532	write(,)"Fixed H2O aerosol distribution?"
533	aerofixh2o=.false. ! default value
534	call getin_p("aerofixh2o",aerofixh2o)
535	write(,)" aerofixh2o = ",aerofixh2o
536
537	write(,)"Radiatively active sulfuric acid aerosols?"
538	aeroh2so4=.false. ! default value
539	call getin_p("aeroh2so4",aeroh2so4)
540	write(,)" aeroh2so4 = ",aeroh2so4
541
542	!=================================
543
544	write(,)"Radiatively active two-layer aerosols?"
545	aeroback2lay=.false. ! default value
546	call getin_p("aeroback2lay",aeroback2lay)
547	write(,)" aeroback2lay = ",aeroback2lay
548
549	write(,)"Radiatively active ammonia cloud?"
550	aeronh3=.false. ! default value
551	call getin_p("aeronh3",aeronh3)
552	write(,)" aeronh3 = ",aeronh3
553
554	write(,)"Radiatively active auroral aerosols?"
555	aeroaurora=.false. ! default value
556	call getin_p("aeroaurora",aeroaurora)
557	write(,)" aeroaurora = ",aeroaurora
558
559	write(,)"TWOLAY AEROSOL: total optical depth ", &
560	"in the tropospheric layer (visible)"
561	obs_tau_col_tropo=8.D0
562	call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
563	write(,)" obs_tau_col_tropo = ",obs_tau_col_tropo
564
565	write(,)"TWOLAY AEROSOL: total optical depth ", &
566	"in the stratospheric layer (visible)"
567	obs_tau_col_strato=0.08D0
568	call getin_p("obs_tau_col_strato",obs_tau_col_strato)
569	write(,)" obs_tau_col_strato = ",obs_tau_col_strato
570
571	write(,)"TWOLAY AEROSOL: optprop_back2lay_vis?"
572	optprop_back2lay_vis = 'optprop_saturn_vis_n20.dat'
573	call getin_p("optprop_back2lay_vis",optprop_back2lay_vis)
574	write(,)" optprop_back2lay_vis = ",optprop_back2lay_vis
575
576	write(,)"TWOLAY AEROSOL: optprop_back2lay_ir?"
577	optprop_back2lay_ir = 'optprop_saturn_ir_n20.dat'
578	call getin_p("optprop_back2lay_ir",optprop_back2lay_ir)
579	write(,)" optprop_back2lay_ir = ",optprop_back2lay_ir
580
581	write(,)"TWOLAY AEROSOL: pres_bottom_tropo? in pa"
582	pres_bottom_tropo=66000.0
583	call getin_p("pres_bottom_tropo",pres_bottom_tropo)
584	write(,)" pres_bottom_tropo = ",pres_bottom_tropo
585
586	write(,)"TWOLAY AEROSOL: pres_top_tropo? in pa"
587	pres_top_tropo=18000.0
588	call getin_p("pres_top_tropo",pres_top_tropo)
589	write(,)" pres_top_tropo = ",pres_top_tropo
590
591	write(,)"TWOLAY AEROSOL: pres_bottom_strato? in pa"
592	pres_bottom_strato=2000.0
593	call getin_p("pres_bottom_strato",pres_bottom_strato)
594	write(,)" pres_bottom_strato = ",pres_bottom_strato
595
596	! Sanity check
597	if (pres_bottom_strato .gt. pres_top_tropo) then
598	print*,'Error : TWOLAY AEROSOL, Please ensure that in callphys.def'
599	print*,'you have pres_top_tropo > pres_bottom_strato !'
600	stop
601	endif
602
603	write(,)"TWOLAY AEROSOL: pres_top_strato? in pa"
604	pres_top_strato=100.0
605	call getin_p("pres_top_strato",pres_top_strato)
606	write(,)" pres_top_strato = ",pres_top_strato
607
608	write(,)"TWOLAY AEROSOL: particle size in the ", &
609	"tropospheric layer, in meters"
610	size_tropo=2.e-6
611	call getin_p("size_tropo",size_tropo)
612	write(,)" size_tropo = ",size_tropo
613
614	write(,)"TWOLAY AEROSOL: particle size in the ", &
615	"stratospheric layer, in meters"
616	size_strato=1.e-7
617	call getin_p("size_strato",size_strato)
618	write(,)" size_strato = ",size_strato
619
620	write(,)"NH3 (thin) cloud: total optical depth"
621	tau_nh3_cloud=7.D0
622	call getin_p("tau_nh3_cloud",tau_nh3_cloud)
623	write(,)" tau_nh3_cloud = ",tau_nh3_cloud
624
625	write(,)"NH3 (thin) cloud pressure level"
626	pres_nh3_cloud=7.D0
627	call getin_p("pres_nh3_cloud",pres_nh3_cloud)
628	write(,)" pres_nh3_cloud = ",pres_nh3_cloud
629
630	write(,)"NH3 (thin) cloud: particle sizes"
631	size_nh3_cloud=3.e-6
632	call getin_p("size_nh3_cloud",size_nh3_cloud)
633	write(,)" size_nh3_cloud = ",size_nh3_cloud
634
635	!=================================
636
637	write(,)"Cloud pressure level (with kastprof only):"
638	cloudlvl=0. ! default value
639	call getin_p("cloudlvl",cloudlvl)
640	write(,)" cloudlvl = ",cloudlvl
641
642	write(,)"Is the variable gas species radiatively active?"
643	Tstrat=167.0
644	varactive=.false.
645	call getin_p("varactive",varactive)
646	write(,)" varactive = ",varactive
647
648	write(,)"Is the variable gas species distribution set?"
649	varfixed=.false.
650	call getin_p("varfixed",varfixed)
651	write(,)" varfixed = ",varfixed
652
653	write(,)"What is the saturation % of the variable species?"
654	satval=0.8
655	call getin_p("satval",satval)
656	write(,)" satval = ",satval
657
658
659	! Test of incompatibility:
660	! if varactive, then varfixed should be false
661	if (varactive.and.varfixed) then
662	print*,'if varactive, varfixed must be OFF!'
663	stop
664	endif
665
666	write(,) "Gravitationnal sedimentation ?"
667	sedimentation=.false. ! default value
668	call getin_p("sedimentation",sedimentation)
669	write(,) " sedimentation = ",sedimentation
670
671	write(,) "Compute water cycle ?"
672	water=.false. ! default value
673	call getin_p("water",water)
674	write(,) " water = ",water
675
676	! Test of incompatibility:
677	! if water is true, there should be at least a tracer
678	if (water.and.(.not.tracer)) then
679	print*,'if water is ON, tracer must be ON too!'
680	stop
681	endif
682
683	write(,) "Include water condensation ?"
684	watercond=.false. ! default value
685	call getin_p("watercond",watercond)
686	write(,) " watercond = ",watercond
687
688	! Test of incompatibility:
689	! if watercond is used, then water should be used too
690	if (watercond.and.(.not.water)) then
691	print*,'if watercond is used, water should be used too'
692	stop
693	endif
694
695	write(,) "Include water precipitation ?"
696	waterrain=.false. ! default value
697	call getin_p("waterrain",waterrain)
698	write(,) " waterrain = ",waterrain
699
700	write(,) "Include surface hydrology ?"
701	hydrology=.false. ! default value
702	call getin_p("hydrology",hydrology)
703	write(,) " hydrology = ",hydrology
704
705	write(,) "Evolve surface water sources ?"
706	sourceevol=.false. ! default value
707	call getin_p("sourceevol",sourceevol)
708	write(,) " sourceevol = ",sourceevol
709
710	write(,) "Ice evolution timestep ?"
711	icetstep=100.0 ! default value
712	call getin_p("icetstep",icetstep)
713	write(,) " icetstep = ",icetstep
714
715	write(,) "Spectral Dependant albedo ?"
716	albedo_spectral_mode=.false. ! default value
717	call getin_p("albedo_spectral_mode",albedo_spectral_mode)
718	write(,) " albedo_spectral_mode = ",albedo_spectral_mode
719
720	write(,) "Snow albedo ?"
721	write(,) "If albedo_spectral_mode=.true., then this "
722	write(,) "corresponds to the 0.5 microns snow albedo."
723	albedosnow=0.5 ! default value
724	call getin_p("albedosnow",albedosnow)
725	write(,) " albedosnow = ",albedosnow
726
727	write(,) "CO2 ice albedo ?"
728	albedoco2ice=0.5 ! default value
729	call getin_p("albedoco2ice",albedoco2ice)
730	write(,) " albedoco2ice = ",albedoco2ice
731
732	write(,) "Maximum ice thickness ?"
733	maxicethick=2.0 ! default value
734	call getin_p("maxicethick",maxicethick)
735	write(,) " maxicethick = ",maxicethick
736
737	write(,) "Freezing point of seawater ?"
738	Tsaldiff=-1.8 ! default value
739	call getin_p("Tsaldiff",Tsaldiff)
740	write(,) " Tsaldiff = ",Tsaldiff
741
742	write(,) "Does user want to force cpp and mugaz?"
743	force_cpp=.false. ! default value
744	call getin_p("force_cpp",force_cpp)
745	write(,) " force_cpp = ",force_cpp
746
747	if (force_cpp) then
748	mugaz = -99999.
749	PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
750	call getin_p("mugaz",mugaz)
751	IF (mugaz.eq.-99999.) THEN
752	PRINT *, "mugaz must be set if force_cpp = T"
753	STOP
754	ELSE
755	write(,) "mugaz=",mugaz
756	ENDIF
757	cpp = -99999.
758	PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
759	call getin_p("cpp",cpp)
760	IF (cpp.eq.-99999.) THEN
761	PRINT *, "cpp must be set if force_cpp = T"
762	STOP
763	ELSE
764	write(,) "cpp=",cpp
765	ENDIF
766	endif ! of if (force_cpp)
767	call su_gases
768	call calc_cpp_mugaz
769
770	PRINT*,'--------------------------------------------'
771	PRINT*
772	PRINT*
773	ELSE
774	write(,)
775	write(,) 'Cannot read file callphys.def. Is it here ?'
776	stop
777	ENDIF ! of IF(iscallphys)
778
779	PRINT*
780	PRINT*,'inifis: daysec',daysec
781	PRINT*
782	PRINT*,'inifis: The radiative transfer is computed:'
783	PRINT*,' each ',iradia,' physical time-step'
784	PRINT,' or each ',iradiadtphys,' seconds'
785	PRINT*
786
787
788	!-----------------------------------------------------------------------
789	! Some more initialization:
790	! ------------------------
791
792	! Initializations for comgeomfi_h
793	#ifndef MESOSCALE
794	totarea=SSUM(ngrid,parea,1)
795	call planetwide_sumval(parea,totarea_planet)
796
797	!! those are defined in comdiurn_h.F90
798	IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
799	IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
800	IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
801	IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
802
803	DO ig=1,ngrid
804	sinlat(ig)=sin(plat(ig))
805	coslat(ig)=cos(plat(ig))
806	sinlon(ig)=sin(plon(ig))
807	coslon(ig)=cos(plon(ig))
808	ENDDO
809	#endif
810	! initialize variables in radinc_h
811	call ini_radinc_h(nlayer)
812
813	! allocate "comsoil_h" arrays
814	call ini_comsoil_h(ngrid)
815
816	END SUBROUTINE inifis
817
818	END MODULE inifis_mod

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

Download in other formats:

Original Format