Context Navigation

inifis_mod.F90 @ 2447

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

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

Download in other formats:

Original Format