Context Navigation

inifis_mod.F90 @ 2299

Last change on this file since 2299 was 2299, checked in by dbardet, 5 years ago

27/04/2017 (r2299) == DB

Add non-orographic gravity waves drag parameterization (Flott scheme adpated from the Earth GCM)

It can be called using 'calllott_nonoro=.true.' in callphys.def, and set the maximum value of the Eliassen-Plam flux 'epflux_max'.
Cumulated output fields are du_nonoro, dv_nonoro (winds tendency due to GW drag), east_gwstress and west_gwstress (stress profile in esatward and westward direction due to GW drag)
These variables are added in Xhisitns and start files.

File size: 31.7 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(,) "Tidally resonant rotation ?"
154	tlocked=.false. ! default value
155	call getin_p("tlocked",tlocked)
156	write(,) "tlocked = ",tlocked
157
158	write(,) "Saturn ring shadowing ?"
159	rings_shadow = .false.
160	call getin_p("rings_shadow", rings_shadow)
161	write(,) "rings_shadow = ", rings_shadow
162
163	write(,) "Compute latitude-dependent gravity field?"
164	oblate = .false.
165	call getin_p("oblate", oblate)
166	write(,) "oblate = ", oblate
167
168	write(,) "Flattening of the planet (a-b)/a "
169	flatten = 0.0
170	call getin_p("flatten", flatten)
171	write(,) "flatten = ", flatten
172
173
174	write(,) "Needed if oblate=.true.: J2"
175	J2 = 0.0
176	call getin_p("J2", J2)
177	write(,) "J2 = ", J2
178
179	write(,) "Needed if oblate=.true.: Planet mass (*1e24 kg)"
180	MassPlanet = 0.0
181	call getin_p("MassPlanet", MassPlanet)
182	write(,) "MassPlanet = ", MassPlanet
183
184	write(,) "Needed if oblate=.true.: Planet mean radius (m)"
185	Rmean = 0.0
186	call getin_p("Rmean", Rmean)
187	write(,) "Rmean = ", Rmean
188
189	! Test of incompatibility:
190	! if tlocked, then diurnal should be false
191	if (tlocked.and.diurnal) then
192	print*,'If diurnal=true, we should turn off tlocked.'
193	stop
194	endif
195
196	write(,) "Tidal resonance ratio ?"
197	nres=0 ! default value
198	call getin_p("nres",nres)
199	write(,) "nres = ",nres
200
201	write(,) "Write some extra output to the screen ?"
202	lwrite=.false. ! default value
203	call getin_p("lwrite",lwrite)
204	write(,) " lwrite = ",lwrite
205
206	write(,) "Save statistics in file stats.nc ?"
207	callstats=.true. ! default value
208	call getin_p("callstats",callstats)
209	write(,) " callstats = ",callstats
210
211	write(,) "Test energy conservation of model physics ?"
212	enertest=.false. ! default value
213	call getin_p("enertest",enertest)
214	write(,) " enertest = ",enertest
215
216	write(,) "Check to see if cpp values used match gases.def ?"
217	check_cpp_match=.true. ! default value
218	call getin_p("check_cpp_match",check_cpp_match)
219	write(,) " check_cpp_match = ",check_cpp_match
220
221	write(,) "call radiative transfer ?"
222	callrad=.true. ! default value
223	call getin_p("callrad",callrad)
224	write(,) " callrad = ",callrad
225
226	write(,) "call correlated-k radiative transfer ?"
227	corrk=.true. ! default value
228	call getin_p("corrk",corrk)
229	write(,) " corrk = ",corrk
230
231	write(,) "prohibit calculations outside corrk T grid?"
232	strictboundcorrk=.true. ! default value
233	call getin_p("strictboundcorrk",strictboundcorrk)
234	write(,) "strictboundcorrk = ",strictboundcorrk
235
236	write(,) "Minimum atmospheric temperature for Planck function integration ?"
237	tplanckmin=30.0 ! default value
238	call getin_p("tplanckmin",tplanckmin)
239	write(,) " tplanckmin = ",tplanckmin
240
241	write(,) "Maximum atmospheric temperature for Planck function integration ?"
242	tplanckmax=1500.0 ! default value
243	call getin_p("tplanckmax",tplanckmax)
244	write(,) " tplanckmax = ",tplanckmax
245
246	write(,) "Temperature step for Planck function integration ?"
247	dtplanck=0.1 ! default value
248	call getin_p("dtplanck",dtplanck)
249	write(,) " dtplanck = ",dtplanck
250
251	write(,) "call gaseous absorption in the visible bands?", &
252	"(matters only if callrad=T)"
253	callgasvis=.false. ! default value
254	call getin_p("callgasvis",callgasvis)
255	write(,) " callgasvis = ",callgasvis
256
257	write(,) "call continuum opacities in radiative transfer ?", &
258	"(matters only if callrad=T)"
259	continuum=.true. ! default value
260	call getin_p("continuum",continuum)
261	write(,) " continuum = ",continuum
262
263	write(,) "use analytic function for H2O continuum ?"
264	H2Ocont_simple=.false. ! default value
265	call getin_p("H2Ocont_simple",H2Ocont_simple)
266	write(,) " H2Ocont_simple = ",H2Ocont_simple
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
430	! Test of incompatibility:
431	! if kastprof used, we must be in 1D
432	if (kastprof.and.(ngrid.gt.1)) then
433	print*,'kastprof can only be used in 1D!'
434	call abort
435	endif
436
437	write(,)"Stratospheric temperature for kastprof mode?"
438	Tstrat=167.0
439	call getin_p("Tstrat",Tstrat)
440	write(,)" Tstrat = ",Tstrat
441
442	write(,)"Remove lower boundary?"
443	nosurf=.false.
444	call getin_p("nosurf",nosurf)
445	write(,)" nosurf = ",nosurf
446
447	! Tests of incompatibility:
448	if (nosurf.and.callsoil) then
449	print*,'nosurf not compatible with soil scheme!'
450	print*,'... got to make a choice!'
451	call abort
452	endif
453
454	write(,)"Add an internal heat flux?", &
455	"... matters only if callsoil=F"
456	intheat=0.
457	call getin_p("intheat",intheat)
458	write(,)" intheat = ",intheat
459
460	write(,)"Use Newtonian cooling for radiative transfer?"
461	newtonian=.false.
462	call getin_p("newtonian",newtonian)
463	write(,)" newtonian = ",newtonian
464
465	! Tests of incompatibility:
466	if (newtonian.and.corrk) then
467	print*,'newtonian not compatible with correlated-k!'
468	call abort
469	endif
470	if (newtonian.and.calladj) then
471	print*,'newtonian not compatible with adjustment!'
472	call abort
473	endif
474	if (newtonian.and.calldifv) then
475	print*,'newtonian not compatible with a boundary layer!'
476	call abort
477	endif
478
479	write(,)"Test physics timescale in 1D?"
480	testradtimes=.false.
481	call getin_p("testradtimes",testradtimes)
482	write(,)" testradtimes = ",testradtimes
483
484	! Test of incompatibility:
485	! if testradtimes used, we must be in 1D
486	if (testradtimes.and.(ngrid.gt.1)) then
487	print*,'testradtimes can only be used in 1D!'
488	call abort
489	endif
490
491	write(,)"Default planetary temperature?"
492	tplanet=215.0
493	call getin_p("tplanet",tplanet)
494	write(,)" tplanet = ",tplanet
495
496	write(,)"Which star?"
497	startype=1 ! default value = Sol
498	call getin_p("startype",startype)
499	write(,)" startype = ",startype
500
501	write(,)"Value of stellar flux at 1 AU?"
502	Fat1AU=1356.0 ! default value = Sol today
503	call getin_p("Fat1AU",Fat1AU)
504	write(,)" Fat1AU = ",Fat1AU
505
506
507	! TRACERS:
508
509	write(,)"Varying H2O cloud fraction?"
510	CLFvarying=.false. ! default value
511	call getin_p("CLFvarying",CLFvarying)
512	write(,)" CLFvarying = ",CLFvarying
513
514	write(,)"Value of fixed H2O cloud fraction?"
515	CLFfixval=1.0 ! default value
516	call getin_p("CLFfixval",CLFfixval)
517	write(,)" CLFfixval = ",CLFfixval
518
519	write(,)"fixed radii for Cloud particles?"
520	radfixed=.false. ! default value
521	call getin_p("radfixed",radfixed)
522	write(,)" radfixed = ",radfixed
523
524	if(kastprof)then
525	radfixed=.true.
526	endif
527
528	write(,)"Number mixing ratio of CO2 ice particles:"
529	Nmix_co2=1.e6 ! default value
530	call getin_p("Nmix_co2",Nmix_co2)
531	write(,)" Nmix_co2 = ",Nmix_co2
532
533	! write(,)"Number of radiatively active aerosols:"
534	! naerkind=0. ! default value
535	! call getin_p("naerkind",naerkind)
536	! write(,)" naerkind = ",naerkind
537
538	write(,)"Opacity of dust (if used):"
539	dusttau=0. ! default value
540	call getin_p("dusttau",dusttau)
541	write(,)" dusttau = ",dusttau
542
543	write(,)"Radiatively active CO2 aerosols?"
544	aeroco2=.false. ! default value
545	call getin_p("aeroco2",aeroco2)
546	write(,)" aeroco2 = ",aeroco2
547
548	write(,)"Fixed CO2 aerosol distribution?"
549	aerofixco2=.false. ! default value
550	call getin_p("aerofixco2",aerofixco2)
551	write(,)" aerofixco2 = ",aerofixco2
552
553	write(,)"Radiatively active water ice?"
554	aeroh2o=.false. ! default value
555	call getin_p("aeroh2o",aeroh2o)
556	write(,)" aeroh2o = ",aeroh2o
557
558	write(,)"Fixed H2O aerosol distribution?"
559	aerofixh2o=.false. ! default value
560	call getin_p("aerofixh2o",aerofixh2o)
561	write(,)" aerofixh2o = ",aerofixh2o
562
563	write(,)"Radiatively active sulfuric acid aerosols?"
564	aeroh2so4=.false. ! default value
565	call getin_p("aeroh2so4",aeroh2so4)
566	write(,)" aeroh2so4 = ",aeroh2so4
567
568	write(,)"Radiatively active auroral aerosols?"
569	aeroaurora=.false. ! default value
570	call getin_p("aeroaurora",aeroaurora)
571	write(,)" aeroaurora = ",aeroaurora
572
573	!=================================
574	! TWOLAY scheme and NH3 cloudare left for retrocompatibility only,
575	! You should now use N-LAYER scheme (see below).
576
577	write(,)"Radiatively active two-layer aerosols?"
578	aeroback2lay=.false. ! default value
579	call getin_p("aeroback2lay",aeroback2lay)
580	write(,)" aeroback2lay = ",aeroback2lay
581
582	if (aeroback2lay) then
583	print*,'Warning : The TWOLAY AEROSOL scheme is deprecated and buggy...'
584	print*,'You should use the generic n-layer scheme (see aeronlay).'
585	endif
586
587	write(,)"Radiatively active ammonia cloud?"
588	aeronh3=.false. ! default value
589	call getin_p("aeronh3",aeronh3)
590	write(,)" aeronh3 = ",aeronh3
591
592	if (aeronh3) then
593	print*,'Warning : You are using specific NH3 cloud scheme ...'
594	print*,'You should use the generic n-layer scheme (see aeronlay).'
595	endif
596
597	write(,)"TWOLAY AEROSOL: total optical depth ", &
598	"in the tropospheric layer (visible)"
599	obs_tau_col_tropo=8.D0
600	call getin_p("obs_tau_col_tropo",obs_tau_col_tropo)
601	write(,)" obs_tau_col_tropo = ",obs_tau_col_tropo
602
603	write(,)"TWOLAY AEROSOL: total optical depth ", &
604	"in the stratospheric layer (visible)"
605	obs_tau_col_strato=0.08D0
606	call getin_p("obs_tau_col_strato",obs_tau_col_strato)
607	write(,)" obs_tau_col_strato = ",obs_tau_col_strato
608
609	write(,)"TWOLAY AEROSOL: optprop_back2lay_vis?"
610	optprop_back2lay_vis = 'optprop_saturn_vis_n20.dat'
611	call getin_p("optprop_back2lay_vis",optprop_back2lay_vis)
612	write(,)" optprop_back2lay_vis = ",optprop_back2lay_vis
613
614	write(,)"TWOLAY AEROSOL: optprop_back2lay_ir?"
615	optprop_back2lay_ir = 'optprop_saturn_ir_n20.dat'
616	call getin_p("optprop_back2lay_ir",optprop_back2lay_ir)
617	write(,)" optprop_back2lay_ir = ",optprop_back2lay_ir
618
619	write(,)"TWOLAY AEROSOL: pres_bottom_tropo? in pa"
620	pres_bottom_tropo=66000.0
621	call getin_p("pres_bottom_tropo",pres_bottom_tropo)
622	write(,)" pres_bottom_tropo = ",pres_bottom_tropo
623
624	write(,)"TWOLAY AEROSOL: pres_top_tropo? in pa"
625	pres_top_tropo=18000.0
626	call getin_p("pres_top_tropo",pres_top_tropo)
627	write(,)" pres_top_tropo = ",pres_top_tropo
628
629	write(,)"TWOLAY AEROSOL: pres_bottom_strato? in pa"
630	pres_bottom_strato=2000.0
631	call getin_p("pres_bottom_strato",pres_bottom_strato)
632	write(,)" pres_bottom_strato = ",pres_bottom_strato
633
634	! Sanity check
635	if (pres_bottom_strato .gt. pres_top_tropo) then
636	print*,'Error : TWOLAY AEROSOL, Please ensure that in callphys.def'
637	print*,'you have pres_top_tropo > pres_bottom_strato !'
638	call abort
639	endif
640
641	write(,)"TWOLAY AEROSOL: pres_top_strato? in pa"
642	pres_top_strato=100.0
643	call getin_p("pres_top_strato",pres_top_strato)
644	write(,)" pres_top_strato = ",pres_top_strato
645
646	write(,)"TWOLAY AEROSOL: particle size in the ", &
647	"tropospheric layer, in meters"
648	size_tropo=2.e-6
649	call getin_p("size_tropo",size_tropo)
650	write(,)" size_tropo = ",size_tropo
651
652	write(,)"TWOLAY AEROSOL: particle size in the ", &
653	"stratospheric layer, in meters"
654	size_strato=1.e-7
655	call getin_p("size_strato",size_strato)
656	write(,)" size_strato = ",size_strato
657
658	write(,)"NH3 (thin) cloud: total optical depth"
659	tau_nh3_cloud=7.D0
660	call getin_p("tau_nh3_cloud",tau_nh3_cloud)
661	write(,)" tau_nh3_cloud = ",tau_nh3_cloud
662
663	write(,)"NH3 (thin) cloud pressure level"
664	pres_nh3_cloud=7.D0
665	call getin_p("pres_nh3_cloud",pres_nh3_cloud)
666	write(,)" pres_nh3_cloud = ",pres_nh3_cloud
667
668	write(,)"NH3 (thin) cloud: particle sizes"
669	size_nh3_cloud=3.e-6
670	call getin_p("size_nh3_cloud",size_nh3_cloud)
671	write(,)" size_nh3_cloud = ",size_nh3_cloud
672
673	!=================================
674	! Generic N-LAYER aerosol scheme
675
676	write(,)"Radiatively active generic n-layer aerosols?"
677	aeronlay=.false. ! default value
678	call getin_p("aeronlay",aeronlay)
679	write(,)" aeronlay = ",aeronlay
680
681	write(,)"Number of generic aerosols layers?"
682	nlayaero=1 ! default value
683	call getin_p("nlayaero",nlayaero)
684	! Avoid to allocate arrays of size 0
685	if (aeronlay .and. nlayaero.lt.1) then
686	print*, " You are trying to set no generic aerosols..."
687	print*, " Set aeronlay=.false. instead ! I abort."
688	call abort
689	endif
690	if (.not. aeronlay) nlayaero=1
691	write(,)" nlayaero = ",nlayaero
692
693	! This is necessary, we just set the number of aerosol layers
694	IF(.NOT.ALLOCATED(aeronlay_tauref)) ALLOCATE(aeronlay_tauref(nlayaero))
695	IF(.NOT.ALLOCATED(aeronlay_lamref)) ALLOCATE(aeronlay_lamref(nlayaero))
696	IF(.NOT.ALLOCATED(aeronlay_choice)) ALLOCATE(aeronlay_choice(nlayaero))
697	IF(.NOT.ALLOCATED(aeronlay_pbot)) ALLOCATE(aeronlay_pbot(nlayaero))
698	IF(.NOT.ALLOCATED(aeronlay_ptop)) ALLOCATE(aeronlay_ptop(nlayaero))
699	IF(.NOT.ALLOCATED(aeronlay_sclhght)) ALLOCATE(aeronlay_sclhght(nlayaero))
700	IF(.NOT.ALLOCATED(aeronlay_size)) ALLOCATE(aeronlay_size(nlayaero))
701	IF(.NOT.ALLOCATED(optprop_aeronlay_ir)) ALLOCATE(optprop_aeronlay_ir(nlayaero))
702	IF(.NOT.ALLOCATED(optprop_aeronlay_vis)) ALLOCATE(optprop_aeronlay_vis(nlayaero))
703
704	write(,)"Generic n-layer aerosols: Optical depth at reference wavelenght"
705	aeronlay_tauref=1.0E-1
706	call getin_p("aeronlay_tauref",aeronlay_tauref)
707	write(,)" aeronlay_tauref = ",aeronlay_tauref
708
709	write(,)"Generic n-layer aerosols: Reference wavelenght for optical depths (m)"
710	aeronlay_lamref=0.6E-6
711	call getin_p("aeronlay_lamref",aeronlay_lamref)
712	write(,)" aeronlay_lamref = ",aeronlay_lamref
713
714	write(,)"Generic n-layer aerosols: Vertical profile choice : &
715	(1) Tau btwn ptop and pbot follows atm. scale height &
716	or (2) Tau above pbot follows its own scale height"
717	aeronlay_choice=1
718	call getin_p("aeronlay_choice",aeronlay_choice)
719	write(,)" aeronlay_choice = ",aeronlay_choice
720
721	write(,)"Generic n-layer aerosols: bottom pressures (Pa)"
722	aeronlay_pbot=2000.0
723	call getin_p("aeronlay_pbot",aeronlay_pbot)
724	write(,)" aeronlay_pbot = ",aeronlay_pbot
725
726	write(,)"Generic n-layer aerosols: (if choice=1) Top pressures (Pa) "
727	aeronlay_ptop=300000.0
728	call getin_p("aeronlay_ptop",aeronlay_ptop)
729	write(,)" aeronlay_ptop = ",aeronlay_ptop
730
731	write(,)"Generic n-layer aerosols: (if choice=2) Scale height / atm. scale height"
732	aeronlay_ptop=0.2
733	call getin_p("aeronlay_sclhght",aeronlay_sclhght)
734	write(,)" aeronlay_sclhght = ",aeronlay_sclhght
735
736	write(,)"Generic n-layer aerosols: particles sizes (m)"
737	aeronlay_size=1.e-6
738	call getin_p("aeronlay_size",aeronlay_size)
739	write(,)" aeronlay_size = ",aeronlay_size
740
741	write(,)"Generic n-layer aerosols: VIS optical properties file"
742	optprop_aeronlay_vis = 'optprop_saturn_vis_n20.dat'
743	call getin_p("optprop_aeronlay_vis",optprop_aeronlay_vis)
744	write(,)" optprop_aeronlay_vis = ",optprop_aeronlay_vis
745
746	write(,)"Generic n-layer aerosols: IR optical properties file"
747	optprop_aeronlay_ir = 'optprop_saturn_ir_n20.dat'
748	call getin_p("optprop_aeronlay_ir",optprop_aeronlay_ir)
749	write(,)" optprop_aeronlay_ir = ",optprop_aeronlay_ir
750
751
752	!=================================
753
754	write(,)"Cloud pressure level (with kastprof only):"
755	cloudlvl=0. ! default value
756	call getin_p("cloudlvl",cloudlvl)
757	write(,)" cloudlvl = ",cloudlvl
758
759	write(,)"Is the variable gas species radiatively active?"
760	Tstrat=167.0
761	varactive=.false.
762	call getin_p("varactive",varactive)
763	write(,)" varactive = ",varactive
764
765	write(,)"Is the variable gas species distribution set?"
766	varfixed=.false.
767	call getin_p("varfixed",varfixed)
768	write(,)" varfixed = ",varfixed
769
770	write(,)"What is the saturation % of the variable species?"
771	satval=0.8
772	call getin_p("satval",satval)
773	write(,)" satval = ",satval
774
775
776	! Test of incompatibility:
777	! if varactive, then varfixed should be false
778	if (varactive.and.varfixed) then
779	print*,'if varactive, varfixed must be OFF!'
780	stop
781	endif
782
783	write(,) "Gravitationnal sedimentation ?"
784	sedimentation=.false. ! default value
785	call getin_p("sedimentation",sedimentation)
786	write(,) " sedimentation = ",sedimentation
787
788	write(,) "Compute water cycle ?"
789	water=.false. ! default value
790	call getin_p("water",water)
791	write(,) " water = ",water
792
793	! Test of incompatibility:
794	! if water is true, there should be at least a tracer
795	if (water.and.(.not.tracer)) then
796	print*,'if water is ON, tracer must be ON too!'
797	stop
798	endif
799
800	write(,) "Include water condensation ?"
801	watercond=.false. ! default value
802	call getin_p("watercond",watercond)
803	write(,) " watercond = ",watercond
804
805	! Test of incompatibility:
806	! if watercond is used, then water should be used too
807	if (watercond.and.(.not.water)) then
808	print*,'if watercond is used, water should be used too'
809	stop
810	endif
811
812	write(,) "Include water precipitation ?"
813	waterrain=.false. ! default value
814	call getin_p("waterrain",waterrain)
815	write(,) " waterrain = ",waterrain
816
817	write(,) "Include surface hydrology ?"
818	hydrology=.false. ! default value
819	call getin_p("hydrology",hydrology)
820	write(,) " hydrology = ",hydrology
821
822	write(,) "Evolve surface water sources ?"
823	sourceevol=.false. ! default value
824	call getin_p("sourceevol",sourceevol)
825	write(,) " sourceevol = ",sourceevol
826
827	write(,) "Ice evolution timestep ?"
828	icetstep=100.0 ! default value
829	call getin_p("icetstep",icetstep)
830	write(,) " icetstep = ",icetstep
831
832	write(,) "Spectral Dependant albedo ?"
833	albedo_spectral_mode=.false. ! default value
834	call getin_p("albedo_spectral_mode",albedo_spectral_mode)
835	write(,) " albedo_spectral_mode = ",albedo_spectral_mode
836
837	write(,) "Snow albedo ?"
838	write(,) "If albedo_spectral_mode=.true., then this "
839	write(,) "corresponds to the 0.5 microns snow albedo."
840	albedosnow=0.5 ! default value
841	call getin_p("albedosnow",albedosnow)
842	write(,) " albedosnow = ",albedosnow
843
844	write(,) "CO2 ice albedo ?"
845	albedoco2ice=0.5 ! default value
846	call getin_p("albedoco2ice",albedoco2ice)
847	write(,) " albedoco2ice = ",albedoco2ice
848
849	write(,) "Maximum ice thickness ?"
850	maxicethick=2.0 ! default value
851	call getin_p("maxicethick",maxicethick)
852	write(,) " maxicethick = ",maxicethick
853
854	write(,) "Freezing point of seawater ?"
855	Tsaldiff=-1.8 ! default value
856	call getin_p("Tsaldiff",Tsaldiff)
857	write(,) " Tsaldiff = ",Tsaldiff
858
859	write(,) "Does user want to force cpp and mugaz?"
860	force_cpp=.false. ! default value
861	call getin_p("force_cpp",force_cpp)
862	write(,) " force_cpp = ",force_cpp
863
864	if (force_cpp) then
865	mugaz = -99999.
866	PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
867	call getin_p("mugaz",mugaz)
868	IF (mugaz.eq.-99999.) THEN
869	PRINT *, "mugaz must be set if force_cpp = T"
870	STOP
871	ELSE
872	write(,) "mugaz=",mugaz
873	ENDIF
874	cpp = -99999.
875	PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
876	call getin_p("cpp",cpp)
877	IF (cpp.eq.-99999.) THEN
878	PRINT *, "cpp must be set if force_cpp = T"
879	STOP
880	ELSE
881	write(,) "cpp=",cpp
882	ENDIF
883	endif ! of if (force_cpp)
884	call su_gases
885	call calc_cpp_mugaz
886
887	PRINT*,'--------------------------------------------'
888	PRINT*
889	PRINT*
890	ELSE
891	write(,)
892	write(,) 'Cannot read file callphys.def. Is it here ?'
893	stop
894	ENDIF ! of IF(iscallphys)
895
896	PRINT*
897	PRINT*,'inifis: daysec',daysec
898	PRINT*
899	PRINT*,'inifis: The radiative transfer is computed:'
900	PRINT*,' each ',iradia,' physical time-step'
901	PRINT,' or each ',iradiadtphys,' seconds'
902	PRINT*
903
904
905	!-----------------------------------------------------------------------
906	! Some more initialization:
907	! ------------------------
908
909	! Initializations for comgeomfi_h
910	#ifndef MESOSCALE
911	totarea=SSUM(ngrid,parea,1)
912	call planetwide_sumval(parea,totarea_planet)
913
914	!! those are defined in comdiurn_h.F90
915	IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
916	IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
917	IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
918	IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
919
920	DO ig=1,ngrid
921	sinlat(ig)=sin(plat(ig))
922	coslat(ig)=cos(plat(ig))
923	sinlon(ig)=sin(plon(ig))
924	coslon(ig)=cos(plon(ig))
925	ENDDO
926	#endif
927	! initialize variables in radinc_h
928	call ini_radinc_h(nlayer,tplanckmin,tplanckmax,dtplanck)
929
930	! allocate "comsoil_h" arrays
931	call ini_comsoil_h(ngrid)
932
933	END SUBROUTINE inifis
934
935	END MODULE inifis_mod

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

Download in other formats:

Original Format