Context Navigation

inifis_mod.F90 @ 2297

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

Add a generic n-layer aerosol scheme to replace the former buggy 2-layer scheme as well as the hard-coded NH3 cloud.

It can be called using 'aeronlay=.true.' in callphys.def, and set the number of layers (up to 4) with 'nlayaero'.
Then, the following parameters are read as arrays of size nlayaero in callphys.def (separated by blank space)

*aeronlay_tauref (Optical depth of aerosol layer at ref wavelenght)
*aeronlay_lamref (Ref wavelenght (m))
*aeronlay_choice (Choice for vertical profile - 1:tau follows atm scale height btwn top and bottom - 2:tau follows it own scale height)
*aeronlay_pbot (Bottom pressure (Pa))
*aeronlay_ptop (Top pressure (Pa) - useful only if choice=1)
*aeronlay_sclhght (Ratio of aerosol layer scale height / atmospheric scale height - useful only if choice=2 )
*aeronlay_size (Particle size (m))
*optprop_aeronlay_vis File for VIS opt properties.
*optprop_aeronlay_ir File for IR opt properties.

+Extra info :

+ In addition of solving the bug from 2-layer it enables different optical properties.
+ The former scheme are left for retrocompatibility (for now) but you should use the new one.
+ See aeropacity.F90 for the calculations

+ Each layer can have different optical properties, size of particle ...
+ You have different choices for vertical profile of the aerosol layers :

aeronlay_choice = 1 : Layer tau is spread between ptop and pbot following atm scale height.
aeronlay_choice = 2 : Layer tau follows its own scale height above cloud deck (pbot).

In this case ptop is dummy and sclhght gives the ratio H_cl/H_atm.

+ The reference wavelenght for input optical depth is now read as input (aeronlay_lamref)
+ Following the last point some comment is added in suaer_corrk about the 'not-really-dummy'ness of IR lamref..

File size: 31.2 KB

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

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

Download in other formats:

Original Format