Context Navigation

inifis_mod.F90 @ 2542

Last change on this file since 2542 was 2542, checked in by aslmd, 3 years ago

Generic GCM:

Large update of the chemical modules

Read chemical network from input files
Init chemistry with ModernTrac?
Photolysis online calculation

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

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

Download in other formats:

Original Format