Context Navigation

inifis_mod.F90 @ 2520

Last change on this file since 2520 was 2520, checked in by Guillaume Chaverot, 4 years ago
update of the H2O continuum
File size: 33.0 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(,)"Production of haze ?"
425	haze=.false. ! default value
426	call getin_p("haze",haze)
427	write(,)" haze = ",haze
428
429	! Global1D mean and solar zenith angle
430
431	if (ngrid.eq.1) then
432	PRINT*, 'Simulate global averaged conditions ?'
433	global1d = .false. ! default value
434	call getin_p("global1d",global1d)
435	write(,) "global1d = ",global1d
436
437	! Test of incompatibility : if global1d is true, there should not be any diurnal cycle.
438	if (global1d.and.diurnal) then
439	call abort_physic("inifis",'if global1d is true, diurnal must be set to false',1)
440	endif
441
442	if (global1d) then
443	PRINT *,'Solar Zenith angle (deg.) ?'
444	PRINT *,'(assumed for averaged solar flux S/4)'
445	szangle=60.0 ! default value
446	call getin_p("szangle",szangle)
447	write(,) "szangle = ",szangle
448	endif
449	endif
450
451	! Test of incompatibility:
452	! if kastprof used, we must be in 1D
453	if (kastprof.and.(ngrid.gt.1)) then
454	print*,'kastprof can only be used in 1D!'
455	call abort
456	endif
457
458	write(,)"Stratospheric temperature for kastprof mode?"
459	Tstrat=167.0
460	call getin_p("Tstrat",Tstrat)
461	write(,)" Tstrat = ",Tstrat
462
463	write(,)"Remove lower boundary?"
464	nosurf=.false.
465	call getin_p("nosurf",nosurf)
466	write(,)" nosurf = ",nosurf
467
468	! Tests of incompatibility:
469	if (nosurf.and.callsoil) then
470	print*,'nosurf not compatible with soil scheme!'
471	print*,'... got to make a choice!'
472	call abort
473	endif
474
475	write(,)"Add an internal heat flux?", &
476	"... matters only if callsoil=F"
477	intheat=0.
478	call getin_p("intheat",intheat)
479	write(,)" intheat = ",intheat
480
481	write(,)"Use Newtonian cooling for radiative transfer?"
482	newtonian=.false.
483	call getin_p("newtonian",newtonian)
484	write(,)" newtonian = ",newtonian
485
486	! Tests of incompatibility:
487	if (newtonian.and.corrk) then
488	print*,'newtonian not compatible with correlated-k!'
489	call abort
490	endif
491	if (newtonian.and.calladj) then
492	print*,'newtonian not compatible with adjustment!'
493	call abort
494	endif
495	if (newtonian.and.calldifv) then
496	print*,'newtonian not compatible with a boundary layer!'
497	call abort
498	endif
499
500	write(,)"Test physics timescale in 1D?"
501	testradtimes=.false.
502	call getin_p("testradtimes",testradtimes)
503	write(,)" testradtimes = ",testradtimes
504
505	! Test of incompatibility:
506	! if testradtimes used, we must be in 1D
507	if (testradtimes.and.(ngrid.gt.1)) then
508	print*,'testradtimes can only be used in 1D!'
509	call abort
510	endif
511
512	write(,)"Default planetary temperature?"
513	tplanet=215.0
514	call getin_p("tplanet",tplanet)
515	write(,)" tplanet = ",tplanet
516
517	write(,)"Which star?"
518	startype=1 ! default value = Sol
519	call getin_p("startype",startype)
520	write(,)" startype = ",startype
521
522	write(,)"Value of stellar flux at 1 AU?"
523	Fat1AU=1356.0 ! default value = Sol today
524	call getin_p("Fat1AU",Fat1AU)
525	write(,)" Fat1AU = ",Fat1AU
526
527
528	! TRACERS:
529
530	write(,)"Varying H2O cloud fraction?"
531	CLFvarying=.false. ! default value
532	call getin_p("CLFvarying",CLFvarying)
533	write(,)" CLFvarying = ",CLFvarying
534
535	write(,)"Value of fixed H2O cloud fraction?"
536	CLFfixval=1.0 ! default value
537	call getin_p("CLFfixval",CLFfixval)
538	write(,)" CLFfixval = ",CLFfixval
539
540	write(,)"fixed radii for Cloud particles?"
541	radfixed=.false. ! default value
542	call getin_p("radfixed",radfixed)
543	write(,)" radfixed = ",radfixed
544
545	if(kastprof)then
546	radfixed=.true.
547	endif
548
549	write(,)"Number mixing ratio of CO2 ice particles:"
550	Nmix_co2=1.e6 ! default value
551	call getin_p("Nmix_co2",Nmix_co2)
552	write(,)" Nmix_co2 = ",Nmix_co2
553
554	! write(,)"Number of radiatively active aerosols:"
555	! naerkind=0. ! default value
556	! call getin_p("naerkind",naerkind)
557	! write(,)" naerkind = ",naerkind
558
559	write(,)"Opacity of dust (if used):"
560	dusttau=0. ! default value
561	call getin_p("dusttau",dusttau)
562	write(,)" dusttau = ",dusttau
563
564	write(,)"Radiatively active CO2 aerosols?"
565	aeroco2=.false. ! default value
566	call getin_p("aeroco2",aeroco2)
567	write(,)" aeroco2 = ",aeroco2
568
569	write(,)"Fixed CO2 aerosol distribution?"
570	aerofixco2=.false. ! default value
571	call getin_p("aerofixco2",aerofixco2)
572	write(,)" aerofixco2 = ",aerofixco2
573
574	write(,)"Radiatively active water ice?"
575	aeroh2o=.false. ! default value
576	call getin_p("aeroh2o",aeroh2o)
577	write(,)" aeroh2o = ",aeroh2o
578
579	write(,)"Fixed H2O aerosol distribution?"
580	aerofixh2o=.false. ! default value
581	call getin_p("aerofixh2o",aerofixh2o)
582	write(,)" aerofixh2o = ",aerofixh2o
583
584	write(,)"Radiatively active sulfuric acid aerosols?"
585	aeroh2so4=.false. ! default value
586	call getin_p("aeroh2so4",aeroh2so4)
587	write(,)" aeroh2so4 = ",aeroh2so4
588
589	write(,)"Radiatively active auroral aerosols?"
590	aeroaurora=.false. ! default value
591	call getin_p("aeroaurora",aeroaurora)
592	write(,)" aeroaurora = ",aeroaurora
593
594	!=================================
595	! TWOLAY scheme and NH3 cloudare left for retrocompatibility only,
596	! You should now use N-LAYER scheme (see below).
597
598	write(,)"Radiatively active two-layer aerosols?"
599	aeroback2lay=.false. ! default value
600	call getin_p("aeroback2lay",aeroback2lay)
601	write(,)" aeroback2lay = ",aeroback2lay
602
603	if (aeroback2lay) then
604	print*,'Warning : The TWOLAY AEROSOL scheme is deprecated and buggy...'
605	print*,'You should use the generic n-layer scheme (see aeronlay).'
606	endif
607
608	write(,)"Radiatively active ammonia cloud?"
609	aeronh3=.false. ! default value
610	call getin_p("aeronh3",aeronh3)
611	write(,)" aeronh3 = ",aeronh3
612
613	if (aeronh3) then
614	print*,'Warning : You are using specific NH3 cloud scheme ...'
615	print*,'You should use the generic n-layer scheme (see aeronlay).'
616	endif
617
618	write(,)"TWOLAY AEROSOL: total optical depth ", &
619	"in the tropospheric layer (visible)"
620	obs_tau_col_tropo=8.D0
621	call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
622	write(,)" obs_tau_col_tropo = ",obs_tau_col_tropo
623
624	write(,)"TWOLAY AEROSOL: total optical depth ", &
625	"in the stratospheric layer (visible)"
626	obs_tau_col_strato=0.08D0
627	call getin_p("obs_tau_col_strato",obs_tau_col_strato)
628	write(,)" obs_tau_col_strato = ",obs_tau_col_strato
629
630	write(,)"TWOLAY AEROSOL: optprop_back2lay_vis?"
631	optprop_back2lay_vis = 'optprop_saturn_vis_n20.dat'
632	call getin_p("optprop_back2lay_vis",optprop_back2lay_vis)
633	write(,)" optprop_back2lay_vis = ",optprop_back2lay_vis
634
635	write(,)"TWOLAY AEROSOL: optprop_back2lay_ir?"
636	optprop_back2lay_ir = 'optprop_saturn_ir_n20.dat'
637	call getin_p("optprop_back2lay_ir",optprop_back2lay_ir)
638	write(,)" optprop_back2lay_ir = ",optprop_back2lay_ir
639
640	write(,)"TWOLAY AEROSOL: pres_bottom_tropo? in pa"
641	pres_bottom_tropo=66000.0
642	call getin_p("pres_bottom_tropo",pres_bottom_tropo)
643	write(,)" pres_bottom_tropo = ",pres_bottom_tropo
644
645	write(,)"TWOLAY AEROSOL: pres_top_tropo? in pa"
646	pres_top_tropo=18000.0
647	call getin_p("pres_top_tropo",pres_top_tropo)
648	write(,)" pres_top_tropo = ",pres_top_tropo
649
650	write(,)"TWOLAY AEROSOL: pres_bottom_strato? in pa"
651	pres_bottom_strato=2000.0
652	call getin_p("pres_bottom_strato",pres_bottom_strato)
653	write(,)" pres_bottom_strato = ",pres_bottom_strato
654
655	! Sanity check
656	if (pres_bottom_strato .gt. pres_top_tropo) then
657	print*,'Error : TWOLAY AEROSOL, Please ensure that in callphys.def'
658	print*,'you have pres_top_tropo > pres_bottom_strato !'
659	call abort
660	endif
661
662	write(,)"TWOLAY AEROSOL: pres_top_strato? in pa"
663	pres_top_strato=100.0
664	call getin_p("pres_top_strato",pres_top_strato)
665	write(,)" pres_top_strato = ",pres_top_strato
666
667	write(,)"TWOLAY AEROSOL: particle size in the ", &
668	"tropospheric layer, in meters"
669	size_tropo=2.e-6
670	call getin_p("size_tropo",size_tropo)
671	write(,)" size_tropo = ",size_tropo
672
673	write(,)"TWOLAY AEROSOL: particle size in the ", &
674	"stratospheric layer, in meters"
675	size_strato=1.e-7
676	call getin_p("size_strato",size_strato)
677	write(,)" size_strato = ",size_strato
678
679	write(,)"NH3 (thin) cloud: total optical depth"
680	tau_nh3_cloud=7.D0
681	call getin_p("tau_nh3_cloud",tau_nh3_cloud)
682	write(,)" tau_nh3_cloud = ",tau_nh3_cloud
683
684	write(,)"NH3 (thin) cloud pressure level"
685	pres_nh3_cloud=7.D0
686	call getin_p("pres_nh3_cloud",pres_nh3_cloud)
687	write(,)" pres_nh3_cloud = ",pres_nh3_cloud
688
689	write(,)"NH3 (thin) cloud: particle sizes"
690	size_nh3_cloud=3.e-6
691	call getin_p("size_nh3_cloud",size_nh3_cloud)
692	write(,)" size_nh3_cloud = ",size_nh3_cloud
693
694	!=================================
695	! Generic N-LAYER aerosol scheme
696
697	write(,)"Radiatively active generic n-layer aerosols?"
698	aeronlay=.false. ! default value
699	call getin_p("aeronlay",aeronlay)
700	write(,)" aeronlay = ",aeronlay
701
702	write(,)"Number of generic aerosols layers?"
703	nlayaero=1 ! default value
704	call getin_p("nlayaero",nlayaero)
705	! Avoid to allocate arrays of size 0
706	if (aeronlay .and. nlayaero.lt.1) then
707	print*, " You are trying to set no generic aerosols..."
708	print*, " Set aeronlay=.false. instead ! I abort."
709	call abort
710	endif
711	if (.not. aeronlay) nlayaero=1
712	write(,)" nlayaero = ",nlayaero
713
714	! This is necessary, we just set the number of aerosol layers
715	IF(.NOT.ALLOCATED(aeronlay_tauref)) ALLOCATE(aeronlay_tauref(nlayaero))
716	IF(.NOT.ALLOCATED(aeronlay_lamref)) ALLOCATE(aeronlay_lamref(nlayaero))
717	IF(.NOT.ALLOCATED(aeronlay_choice)) ALLOCATE(aeronlay_choice(nlayaero))
718	IF(.NOT.ALLOCATED(aeronlay_pbot)) ALLOCATE(aeronlay_pbot(nlayaero))
719	IF(.NOT.ALLOCATED(aeronlay_ptop)) ALLOCATE(aeronlay_ptop(nlayaero))
720	IF(.NOT.ALLOCATED(aeronlay_sclhght)) ALLOCATE(aeronlay_sclhght(nlayaero))
721	IF(.NOT.ALLOCATED(aeronlay_size)) ALLOCATE(aeronlay_size(nlayaero))
722	IF(.NOT.ALLOCATED(aeronlay_nueff)) ALLOCATE(aeronlay_nueff(nlayaero))
723	IF(.NOT.ALLOCATED(optprop_aeronlay_ir)) ALLOCATE(optprop_aeronlay_ir(nlayaero))
724	IF(.NOT.ALLOCATED(optprop_aeronlay_vis)) ALLOCATE(optprop_aeronlay_vis(nlayaero))
725
726	write(,)"Generic n-layer aerosols: Optical depth at reference wavelenght"
727	aeronlay_tauref=1.0E-1
728	call getin_p("aeronlay_tauref",aeronlay_tauref)
729	write(,)" aeronlay_tauref = ",aeronlay_tauref
730
731	write(,)"Generic n-layer aerosols: Reference wavelenght for optical depths (m)"
732	aeronlay_lamref=0.6E-6
733	call getin_p("aeronlay_lamref",aeronlay_lamref)
734	write(,)" aeronlay_lamref = ",aeronlay_lamref
735
736	write(,)"Generic n-layer aerosols: Vertical profile choice : &
737	(1) Tau btwn ptop and pbot follows atm. scale height &
738	or (2) Tau above pbot follows its own scale height"
739	aeronlay_choice=1
740	call getin_p("aeronlay_choice",aeronlay_choice)
741	write(,)" aeronlay_choice = ",aeronlay_choice
742
743	write(,)"Generic n-layer aerosols: bottom pressures (Pa)"
744	aeronlay_pbot=2000.0
745	call getin_p("aeronlay_pbot",aeronlay_pbot)
746	write(,)" aeronlay_pbot = ",aeronlay_pbot
747
748	write(,)"Generic n-layer aerosols: (if choice=1) Top pressures (Pa) "
749	aeronlay_ptop=300000.0
750	call getin_p("aeronlay_ptop",aeronlay_ptop)
751	write(,)" aeronlay_ptop = ",aeronlay_ptop
752
753	write(,)"Generic n-layer aerosols: (if choice=2) Scale height / atm. scale height"
754	aeronlay_sclhght=0.2
755	call getin_p("aeronlay_sclhght",aeronlay_sclhght)
756	write(,)" aeronlay_sclhght = ",aeronlay_sclhght
757
758	write(,)"Generic n-layer aerosols: particles effective radii (m)"
759	aeronlay_size=1.e-6
760	call getin_p("aeronlay_size",aeronlay_size)
761	write(,)" aeronlay_size = ",aeronlay_size
762
763	write(,)"Generic n-layer aerosols: particles radii effective variance"
764	aeronlay_nueff=0.1
765	call getin_p("aeronlay_nueff",aeronlay_nueff)
766	write(,)" aeronlay_nueff = ",aeronlay_nueff
767
768	write(,)"Generic n-layer aerosols: VIS optical properties file"
769	optprop_aeronlay_vis = 'optprop_saturn_vis_n20.dat'
770	call getin_p("optprop_aeronlay_vis",optprop_aeronlay_vis)
771	write(,)" optprop_aeronlay_vis = ",optprop_aeronlay_vis
772
773	write(,)"Generic n-layer aerosols: IR optical properties file"
774	optprop_aeronlay_ir = 'optprop_saturn_ir_n20.dat'
775	call getin_p("optprop_aeronlay_ir",optprop_aeronlay_ir)
776	write(,)" optprop_aeronlay_ir = ",optprop_aeronlay_ir
777
778
779	!=================================
780
781	write(,)"Cloud pressure level (with kastprof only):"
782	cloudlvl=0. ! default value
783	call getin_p("cloudlvl",cloudlvl)
784	write(,)" cloudlvl = ",cloudlvl
785
786	write(,)"Is the variable gas species radiatively active?"
787	Tstrat=167.0
788	varactive=.false.
789	call getin_p("varactive",varactive)
790	write(,)" varactive = ",varactive
791
792	write(,)"Is the variable gas species distribution set?"
793	varfixed=.false.
794	call getin_p("varfixed",varfixed)
795	write(,)" varfixed = ",varfixed
796
797	write(,)"What is the saturation % of the variable species?"
798	satval=0.8
799	call getin_p("satval",satval)
800	write(,)" satval = ",satval
801
802
803	! Test of incompatibility:
804	! if varactive, then varfixed should be false
805	if (varactive.and.varfixed) then
806	print*,'if varactive, varfixed must be OFF!'
807	stop
808	endif
809
810	write(,) "Gravitationnal sedimentation ?"
811	sedimentation=.false. ! default value
812	call getin_p("sedimentation",sedimentation)
813	write(,) " sedimentation = ",sedimentation
814
815	write(,) "Compute water cycle ?"
816	water=.false. ! default value
817	call getin_p("water",water)
818	write(,) " water = ",water
819
820	! Test of incompatibility:
821	! if water is true, there should be at least a tracer
822	if (water.and.(.not.tracer)) then
823	print*,'if water is ON, tracer must be ON too!'
824	stop
825	endif
826
827	write(,) "Include water condensation ?"
828	watercond=.false. ! default value
829	call getin_p("watercond",watercond)
830	write(,) " watercond = ",watercond
831
832	! Test of incompatibility:
833	! if watercond is used, then water should be used too
834	if (watercond.and.(.not.water)) then
835	print*,'if watercond is used, water should be used too'
836	stop
837	endif
838
839	write(,) "Include water precipitation ?"
840	waterrain=.false. ! default value
841	call getin_p("waterrain",waterrain)
842	write(,) " waterrain = ",waterrain
843
844	write(,) "Include surface hydrology ?"
845	hydrology=.false. ! default value
846	call getin_p("hydrology",hydrology)
847	write(,) " hydrology = ",hydrology
848
849	write(,) "Evolve surface water sources ?"
850	sourceevol=.false. ! default value
851	call getin_p("sourceevol",sourceevol)
852	write(,) " sourceevol = ",sourceevol
853
854	write(,) "Ice evolution timestep ?"
855	icetstep=100.0 ! default value
856	call getin_p("icetstep",icetstep)
857	write(,) " icetstep = ",icetstep
858
859	write(,) "Spectral Dependant albedo ?"
860	albedo_spectral_mode=.false. ! default value
861	call getin_p("albedo_spectral_mode",albedo_spectral_mode)
862	write(,) " albedo_spectral_mode = ",albedo_spectral_mode
863
864	write(,) "Snow albedo ?"
865	write(,) "If albedo_spectral_mode=.true., then this "
866	write(,) "corresponds to the 0.5 microns snow albedo."
867	albedosnow=0.5 ! default value
868	call getin_p("albedosnow",albedosnow)
869	write(,) " albedosnow = ",albedosnow
870
871	write(,) "Ocean albedo ?"
872	alb_ocean=0.07 ! default value
873	call getin_p("alb_ocean",alb_ocean)
874	write(,) " alb_ocean = ",alb_ocean
875
876	write(,) "CO2 ice albedo ?"
877	albedoco2ice=0.5 ! default value
878	call getin_p("albedoco2ice",albedoco2ice)
879	write(,) " albedoco2ice = ",albedoco2ice
880
881	write(,) "Maximum ice thickness ?"
882	maxicethick=2.0 ! default value
883	call getin_p("maxicethick",maxicethick)
884	write(,) " maxicethick = ",maxicethick
885
886	write(,) "Freezing point of seawater ?"
887	Tsaldiff=-1.8 ! default value
888	call getin_p("Tsaldiff",Tsaldiff)
889	write(,) " Tsaldiff = ",Tsaldiff
890
891	write(,) "Minimum eddy mix coeff in 1D ?"
892	kmixmin=1.0e-2 ! default value
893	call getin_p("kmixmin",kmixmin)
894	write(,) " kmixmin = ",kmixmin
895
896	write(,) "Does user want to force cpp and mugaz?"
897	force_cpp=.false. ! default value
898	call getin_p("force_cpp",force_cpp)
899	write(,) " force_cpp = ",force_cpp
900
901	if (force_cpp) then
902	mugaz = -99999.
903	PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
904	call getin_p("mugaz",mugaz)
905	IF (mugaz.eq.-99999.) THEN
906	PRINT *, "mugaz must be set if force_cpp = T"
907	STOP
908	ELSE
909	write(,) "mugaz=",mugaz
910	ENDIF
911	cpp = -99999.
912	PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
913	call getin_p("cpp",cpp)
914	IF (cpp.eq.-99999.) THEN
915	PRINT *, "cpp must be set if force_cpp = T"
916	STOP
917	ELSE
918	write(,) "cpp=",cpp
919	ENDIF
920	endif ! of if (force_cpp)
921	call su_gases
922	call calc_cpp_mugaz
923
924	PRINT*,'--------------------------------------------'
925	PRINT*
926	PRINT*
927	ELSE
928	write(,)
929	write(,) 'Cannot read file callphys.def. Is it here ?'
930	stop
931	ENDIF ! of IF(iscallphys)
932
933	PRINT*
934	PRINT*,'inifis: daysec',daysec
935	PRINT*
936	PRINT*,'inifis: The radiative transfer is computed:'
937	PRINT*,' each ',iradia,' physical time-step'
938	PRINT,' or each ',iradiadtphys,' seconds'
939	PRINT*
940
941
942	!-----------------------------------------------------------------------
943	! Some more initialization:
944	! ------------------------
945
946	! Initializations for comgeomfi_h
947	#ifndef MESOSCALE
948	totarea=SSUM(ngrid,parea,1)
949	call planetwide_sumval(parea,totarea_planet)
950
951	!! those are defined in comdiurn_h.F90
952	IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
953	IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
954	IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
955	IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
956
957	DO ig=1,ngrid
958	sinlat(ig)=sin(plat(ig))
959	coslat(ig)=cos(plat(ig))
960	sinlon(ig)=sin(plon(ig))
961	coslon(ig)=cos(plon(ig))
962	ENDDO
963	#endif
964	! initialize variables in radinc_h
965	call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
966
967	! allocate "comsoil_h" arrays
968	call ini_comsoil_h(ngrid)
969
970	END SUBROUTINE inifis
971
972	END MODULE inifis_mod

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

Download in other formats:

Original Format