Context Navigation

inifis_mod.F90 @ 2309

Last change on this file since 2309 was 2309, checked in by jvatant, 5 years ago
Fix a nasty copy-paste bug from r2297 in n-layer aerosol scheme --JVO
File size: 31.9 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(optprop_aeronlay_ir)) ALLOCATE(optprop_aeronlay_ir(nlayaero))
707	IF(.NOT.ALLOCATED(optprop_aeronlay_vis)) ALLOCATE(optprop_aeronlay_vis(nlayaero))
708
709	write(,)"Generic n-layer aerosols: Optical depth at reference wavelenght"
710	aeronlay_tauref=1.0E-1
711	call getin_p("aeronlay_tauref",aeronlay_tauref)
712	write(,)" aeronlay_tauref = ",aeronlay_tauref
713
714	write(,)"Generic n-layer aerosols: Reference wavelenght for optical depths (m)"
715	aeronlay_lamref=0.6E-6
716	call getin_p("aeronlay_lamref",aeronlay_lamref)
717	write(,)" aeronlay_lamref = ",aeronlay_lamref
718
719	write(,)"Generic n-layer aerosols: Vertical profile choice : &
720	(1) Tau btwn ptop and pbot follows atm. scale height &
721	or (2) Tau above pbot follows its own scale height"
722	aeronlay_choice=1
723	call getin_p("aeronlay_choice",aeronlay_choice)
724	write(,)" aeronlay_choice = ",aeronlay_choice
725
726	write(,)"Generic n-layer aerosols: bottom pressures (Pa)"
727	aeronlay_pbot=2000.0
728	call getin_p("aeronlay_pbot",aeronlay_pbot)
729	write(,)" aeronlay_pbot = ",aeronlay_pbot
730
731	write(,)"Generic n-layer aerosols: (if choice=1) Top pressures (Pa) "
732	aeronlay_ptop=300000.0
733	call getin_p("aeronlay_ptop",aeronlay_ptop)
734	write(,)" aeronlay_ptop = ",aeronlay_ptop
735
736	write(,)"Generic n-layer aerosols: (if choice=2) Scale height / atm. scale height"
737	aeronlay_sclhght=0.2
738	call getin_p("aeronlay_sclhght",aeronlay_sclhght)
739	write(,)" aeronlay_sclhght = ",aeronlay_sclhght
740
741	write(,)"Generic n-layer aerosols: particles sizes (m)"
742	aeronlay_size=1.e-6
743	call getin_p("aeronlay_size",aeronlay_size)
744	write(,)" aeronlay_size = ",aeronlay_size
745
746	write(,)"Generic n-layer aerosols: VIS optical properties file"
747	optprop_aeronlay_vis = 'optprop_saturn_vis_n20.dat'
748	call getin_p("optprop_aeronlay_vis",optprop_aeronlay_vis)
749	write(,)" optprop_aeronlay_vis = ",optprop_aeronlay_vis
750
751	write(,)"Generic n-layer aerosols: IR optical properties file"
752	optprop_aeronlay_ir = 'optprop_saturn_ir_n20.dat'
753	call getin_p("optprop_aeronlay_ir",optprop_aeronlay_ir)
754	write(,)" optprop_aeronlay_ir = ",optprop_aeronlay_ir
755
756
757	!=================================
758
759	write(,)"Cloud pressure level (with kastprof only):"
760	cloudlvl=0. ! default value
761	call getin_p("cloudlvl",cloudlvl)
762	write(,)" cloudlvl = ",cloudlvl
763
764	write(,)"Is the variable gas species radiatively active?"
765	Tstrat=167.0
766	varactive=.false.
767	call getin_p("varactive",varactive)
768	write(,)" varactive = ",varactive
769
770	write(,)"Is the variable gas species distribution set?"
771	varfixed=.false.
772	call getin_p("varfixed",varfixed)
773	write(,)" varfixed = ",varfixed
774
775	write(,)"What is the saturation % of the variable species?"
776	satval=0.8
777	call getin_p("satval",satval)
778	write(,)" satval = ",satval
779
780
781	! Test of incompatibility:
782	! if varactive, then varfixed should be false
783	if (varactive.and.varfixed) then
784	print*,'if varactive, varfixed must be OFF!'
785	stop
786	endif
787
788	write(,) "Gravitationnal sedimentation ?"
789	sedimentation=.false. ! default value
790	call getin_p("sedimentation",sedimentation)
791	write(,) " sedimentation = ",sedimentation
792
793	write(,) "Compute water cycle ?"
794	water=.false. ! default value
795	call getin_p("water",water)
796	write(,) " water = ",water
797
798	! Test of incompatibility:
799	! if water is true, there should be at least a tracer
800	if (water.and.(.not.tracer)) then
801	print*,'if water is ON, tracer must be ON too!'
802	stop
803	endif
804
805	write(,) "Include water condensation ?"
806	watercond=.false. ! default value
807	call getin_p("watercond",watercond)
808	write(,) " watercond = ",watercond
809
810	! Test of incompatibility:
811	! if watercond is used, then water should be used too
812	if (watercond.and.(.not.water)) then
813	print*,'if watercond is used, water should be used too'
814	stop
815	endif
816
817	write(,) "Include water precipitation ?"
818	waterrain=.false. ! default value
819	call getin_p("waterrain",waterrain)
820	write(,) " waterrain = ",waterrain
821
822	write(,) "Include surface hydrology ?"
823	hydrology=.false. ! default value
824	call getin_p("hydrology",hydrology)
825	write(,) " hydrology = ",hydrology
826
827	write(,) "Evolve surface water sources ?"
828	sourceevol=.false. ! default value
829	call getin_p("sourceevol",sourceevol)
830	write(,) " sourceevol = ",sourceevol
831
832	write(,) "Ice evolution timestep ?"
833	icetstep=100.0 ! default value
834	call getin_p("icetstep",icetstep)
835	write(,) " icetstep = ",icetstep
836
837	write(,) "Spectral Dependant albedo ?"
838	albedo_spectral_mode=.false. ! default value
839	call getin_p("albedo_spectral_mode",albedo_spectral_mode)
840	write(,) " albedo_spectral_mode = ",albedo_spectral_mode
841
842	write(,) "Snow albedo ?"
843	write(,) "If albedo_spectral_mode=.true., then this "
844	write(,) "corresponds to the 0.5 microns snow albedo."
845	albedosnow=0.5 ! default value
846	call getin_p("albedosnow",albedosnow)
847	write(,) " albedosnow = ",albedosnow
848
849	write(,) "CO2 ice albedo ?"
850	albedoco2ice=0.5 ! default value
851	call getin_p("albedoco2ice",albedoco2ice)
852	write(,) " albedoco2ice = ",albedoco2ice
853
854	write(,) "Maximum ice thickness ?"
855	maxicethick=2.0 ! default value
856	call getin_p("maxicethick",maxicethick)
857	write(,) " maxicethick = ",maxicethick
858
859	write(,) "Freezing point of seawater ?"
860	Tsaldiff=-1.8 ! default value
861	call getin_p("Tsaldiff",Tsaldiff)
862	write(,) " Tsaldiff = ",Tsaldiff
863
864	write(,) "Does user want to force cpp and mugaz?"
865	force_cpp=.false. ! default value
866	call getin_p("force_cpp",force_cpp)
867	write(,) " force_cpp = ",force_cpp
868
869	if (force_cpp) then
870	mugaz = -99999.
871	PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
872	call getin_p("mugaz",mugaz)
873	IF (mugaz.eq.-99999.) THEN
874	PRINT *, "mugaz must be set if force_cpp = T"
875	STOP
876	ELSE
877	write(,) "mugaz=",mugaz
878	ENDIF
879	cpp = -99999.
880	PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
881	call getin_p("cpp",cpp)
882	IF (cpp.eq.-99999.) THEN
883	PRINT *, "cpp must be set if force_cpp = T"
884	STOP
885	ELSE
886	write(,) "cpp=",cpp
887	ENDIF
888	endif ! of if (force_cpp)
889	call su_gases
890	call calc_cpp_mugaz
891
892	PRINT*,'--------------------------------------------'
893	PRINT*
894	PRINT*
895	ELSE
896	write(,)
897	write(,) 'Cannot read file callphys.def. Is it here ?'
898	stop
899	ENDIF ! of IF(iscallphys)
900
901	PRINT*
902	PRINT*,'inifis: daysec',daysec
903	PRINT*
904	PRINT*,'inifis: The radiative transfer is computed:'
905	PRINT*,' each ',iradia,' physical time-step'
906	PRINT,' or each ',iradiadtphys,' seconds'
907	PRINT*
908
909
910	!-----------------------------------------------------------------------
911	! Some more initialization:
912	! ------------------------
913
914	! Initializations for comgeomfi_h
915	#ifndef MESOSCALE
916	totarea=SSUM(ngrid,parea,1)
917	call planetwide_sumval(parea,totarea_planet)
918
919	!! those are defined in comdiurn_h.F90
920	IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
921	IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
922	IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
923	IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
924
925	DO ig=1,ngrid
926	sinlat(ig)=sin(plat(ig))
927	coslat(ig)=cos(plat(ig))
928	sinlon(ig)=sin(plon(ig))
929	coslon(ig)=cos(plon(ig))
930	ENDDO
931	#endif
932	! initialize variables in radinc_h
933	call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
934
935	! allocate "comsoil_h" arrays
936	call ini_comsoil_h(ngrid)
937
938	END SUBROUTINE inifis
939
940	END MODULE inifis_mod

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

Download in other formats:

Original Format