Context Navigation

conf_phys.F @ 2162

Last change on this file since 2162 was 2160, checked in by mvals, 6 years ago

Mars GCM:
Set adpatable parameters for the rocket dust storm scheme (parameters included in callkeys.h, and adaptable according to the callphys.def with the function "call getin" in conf_phys.F):

ti_injection, tf_injection (by default: ti_injection=10. and tf_injection=12., impacted files: compute_dtau_mod.F90, vdifc_mod.F)
coeff_detrainment (by default: coeff_detrainment=0., impacted files: rocketduststorm_mod.F90)

File size: 33.4 KB

Line
1	SUBROUTINE conf_phys(ngrid,nlayer,nq)
2
3	!=======================================================================
4	!
5	! purpose:
6	! -------
7	!
8	! Initialisation for the physical parametrisations of the LMD
9	! martian atmospheric general circulation modele.
10	!
11	! author: Frederic Hourdin 15 / 10 /93
12	! -------
13	! modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
14	! Ehouarn Millour (oct. 2008) tracers are now identified
15	! by their names and may not be contiguously
16	! stored in the q(:,:,:,:) array
17	! E.M. (june 2009) use getin routine to load parameters
18	! adapted to the mesoscale use - Aymeric Spiga - 01/2007-07/2011
19	! separated inifis into conf_phys and phys_state_var_init (A. Spiga)
20	!
21	!
22	! arguments:
23	! ----------
24	!
25	! input:
26	! ------
27	!
28	! nq Number of tracers
29	!
30	!=======================================================================
31	!
32	!-----------------------------------------------------------------------
33	! declarations:
34	! -------------
35	! to use 'getin'
36	USE ioipsl_getincom, only : getin
37	use tracer_mod, only : nuice_sed, ccn_factor, nuiceco2_sed,
38	& nuice_ref,nuiceco2_ref
39	use surfdat_h, only: albedo_h2o_ice, inert_h2o_ice,
40	& frost_albedo_threshold
41	use time_phylmdz_mod, only: ecritphy,day_step,iphysiq,ecritstart,
42	& daysec,dtphys
43	use dimradmars_mod, only: naerkind, name_iaer,
44	& ini_scatterers,tauvis
45	use datafile_mod, only: datadir
46
47	IMPLICIT NONE
48	include "callkeys.h"
49	include "microphys.h"
50
51	INTEGER,INTENT(IN) :: ngrid,nlayer,nq
52	INTEGER ierr,j
53
54	CHARACTER ch1*12
55	#ifndef MESOSCALE
56	! read in some parameters from "run.def" for physics,
57	! or shared between dynamics and physics.
58	ecritphy=240 ! default value
59	call getin("ecritphy",ecritphy) ! frequency of outputs in physics,
60	! in dynamical steps
61	day_step=960 ! default value
62	call getin("day_step",day_step) ! number of dynamical steps per day
63	iphysiq=20 ! default value
64	call getin("iphysiq",iphysiq) ! call physics every iphysiq dyn step
65	ecritstart=0 ! default value
66	call getin("ecritstart",ecritstart) ! write a restart every ecristart steps
67	#endif
68
69	! --------------------------------------------------------------
70	! Reading the "callphys.def" file controlling some key options
71	! --------------------------------------------------------------
72
73	! check that 'callphys.def' file is around
74	OPEN(99,file='callphys.def',status='old',form='formatted'
75	& ,iostat=ierr)
76	CLOSE(99)
77
78	IF(ierr.EQ.0) THEN
79	PRINT*
80	PRINT*
81	PRINT*,'--------------------------------------------'
82	PRINT*,' conf_phys: Parameters for the physics (callphys.def)'
83	PRINT*,'--------------------------------------------'
84
85	write(,) "Directory where external input files are:"
86	! default path is set in datafile_mod
87	call getin("datadir",datadir)
88	write(,) " datadir = ",trim(datadir)
89
90	write(,) "Run with or without tracer transport ?"
91	tracer=.false. ! default value
92	call getin("tracer",tracer)
93	write(,) " tracer = ",tracer
94
95	write(,) "Diurnal cycle ?"
96	write(,) "(if diurnal=False, diurnal averaged solar heating)"
97	diurnal=.true. ! default value
98	call getin("diurnal",diurnal)
99	write(,) " diurnal = ",diurnal
100
101	write(,) "Seasonal cycle ?"
102	write(,) "(if season=False, Ls stays constant, to value ",
103	& "set in 'start'"
104	season=.true. ! default value
105	call getin("season",season)
106	write(,) " season = ",season
107
108	write(,) "Write some extra output to the screen ?"
109	lwrite=.false. ! default value
110	call getin("lwrite",lwrite)
111	write(,) " lwrite = ",lwrite
112
113	write(,) "Save statistics in file stats.nc ?"
114	#ifdef MESOSCALE
115	callstats=.false. ! default value
116	#else
117	callstats=.true. ! default value
118	#endif
119	call getin("callstats",callstats)
120	write(,) " callstats = ",callstats
121
122	write(,) "Save EOF profiles in file 'profiles' for ",
123	& "Climate Database?"
124	calleofdump=.false. ! default value
125	call getin("calleofdump",calleofdump)
126	write(,) " calleofdump = ",calleofdump
127
128	write(,) "Dust scenario: 1=constant dust (read from startfi",
129	& " or set as tauvis); 2=Viking scenario; =3 MGS scenario,",
130	& "=6 cold (low dust) scenario; =7 warm (high dust) scenario ",
131	& "=24,25 ... 30 :Mars Year 24, ... or 30 from TES assimilation"
132	iaervar=3 ! default value
133	call getin("iaervar",iaervar)
134	write(,) " iaervar = ",iaervar
135
136	write(,) "Reference (visible) dust opacity at 610 Pa ",
137	& "(matters only if iaervar=1)"
138	! NB: default value of tauvis is set/read in startfi.nc file
139	call getin("tauvis",tauvis)
140	write(,) " tauvis = ",tauvis
141
142	write(,) "Dust vertical distribution:"
143	write(,) "(=1 top set by topdustref parameter;",
144	& " =2 Viking scenario; =3 MGS scenario)"
145	iddist=3 ! default value
146	call getin("iddist",iddist)
147	write(,) " iddist = ",iddist
148
149	write(,) "Dust top altitude (km). (Matters only if iddist=1)"
150	topdustref= 90.0 ! default value
151	call getin("topdustref",topdustref)
152	write(,) " topdustref = ",topdustref
153
154	write(,) "Prescribed surface thermal flux (H/(rho*cp),K m/s)"
155	tke_heat_flux=0. ! default value
156	call getin("tke_heat_flux",tke_heat_flux)
157	write(,) " tke_heat_flux = ",tke_heat_flux
158	write(,) " 0 means the usual schemes are computing"
159
160	write(,) "call radiative transfer ?"
161	callrad=.true. ! default value
162	call getin("callrad",callrad)
163	write(,) " callrad = ",callrad
164
165	write(,) "call slope insolation scheme ?",
166	& "(matters only if callrad=T)"
167	#ifdef MESOSCALE
168	callslope=.true. ! default value
169	#else
170	callslope=.false. ! default value (not supported yet)
171	#endif
172	call getin("callslope",callslope)
173	write(,) " callslope = ",callslope
174
175	write(,) "call NLTE radiative schemes ?",
176	& "(matters only if callrad=T)"
177	callnlte=.false. ! default value
178	call getin("callnlte",callnlte)
179	write(,) " callnlte = ",callnlte
180
181	nltemodel=0 !default value
182	write(,) "NLTE model?"
183	write(,) "0 -> old model, static O"
184	write(,) "1 -> old model, dynamic O"
185	write(,) "2 -> new model"
186	write(,) "(matters only if callnlte=T)"
187	call getin("nltemodel",nltemodel)
188	write(,) " nltemodel = ",nltemodel
189
190	write(,) "call CO2 NIR absorption ?",
191	& "(matters only if callrad=T)"
192	callnirco2=.false. ! default value
193	call getin("callnirco2",callnirco2)
194	write(,) " callnirco2 = ",callnirco2
195
196	write(,) "New NIR NLTE correction ?",
197	$ "0-> old model (no correction)",
198	$ "1-> new correction",
199	$ "(matters only if callnirco2=T)"
200	#ifdef MESOSCALE
201	nircorr=0 !default value. this is OK below 60 km.
202	#else
203	nircorr=0 !default value
204	#endif
205	call getin("nircorr",nircorr)
206	write(,) " nircorr = ",nircorr
207
208	write(,) "call turbulent vertical diffusion ?"
209	calldifv=.true. ! default value
210	call getin("calldifv",calldifv)
211	write(,) " calldifv = ",calldifv
212
213	write(,) "call thermals ?"
214	calltherm=.false. ! default value
215	call getin("calltherm",calltherm)
216	write(,) " calltherm = ",calltherm
217
218	write(,) "call convective adjustment ?"
219	calladj=.true. ! default value
220	call getin("calladj",calladj)
221	write(,) " calladj = ",calladj
222
223	if (calltherm .and. calladj) then
224	print*,'!!! PLEASE NOTE !!!'
225	print*,'convective adjustment is on'
226	print*,'but since thermal plume model is on'
227	print*,'convadj is only activated above the PBL'
228	endif
229
230	write(,) "used latest version of yamada scheme?"
231	callyamada4=.true. ! default value
232	call getin("callyamada4",callyamada4)
233	write(,) " callyamada4 = ",callyamada4
234
235	if (calltherm .and. .not.callyamada4) then
236	print*,'!!!! WARNING WARNING WARNING !!!!'
237	print*,'if calltherm=T we strongly advise that '
238	print*,'you set the flag callyamada4 to T '
239	print*,'!!!! WARNING WARNING WARNING !!!!'
240	endif
241
242	write(,) "call Richardson-based surface layer ?"
243	callrichsl=.false. ! default value
244	call getin("callrichsl",callrichsl)
245	write(,) " callrichsl = ",callrichsl
246
247	if (calltherm .and. .not.callrichsl) then
248	print*,'WARNING WARNING WARNING'
249	print*,'if calltherm=T we strongly advise that '
250	print*,'you use the new surface layer scheme '
251	print*,'by setting callrichsl=T '
252	endif
253
254	if (calladj .and. callrichsl .and. (.not. calltherm)) then
255	print*,'You should not be calling the convective adjustment
256	& scheme with the Richardson surface-layer and without the thermals
257	&. This approach is not
258	& physically consistent and can lead to unrealistic friction
259	& values.'
260	print*,'If you want to use the Ri. surface-layer, either
261	& activate thermals OR de-activate the convective adjustment.'
262	stop
263	endif
264
265	write(,) "call CO2 condensation ?"
266	callcond=.true. ! default value
267	call getin("callcond",callcond)
268	write(,) " callcond = ",callcond
269
270	write(,)"call thermal conduction in the soil ?"
271	callsoil=.true. ! default value
272	call getin("callsoil",callsoil)
273	write(,) " callsoil = ",callsoil
274
275
276	write(,)"call Lott's gravity wave/subgrid topography ",
277	& "scheme ?"
278	calllott=.true. ! default value
279	call getin("calllott",calllott)
280	write(,)" calllott = ",calllott
281
282	write(,)"call Lott's non-oro GWs parameterisation ",
283	& "scheme ?"
284	calllott_nonoro=.false. ! default value
285	call getin("calllott_nonoro",calllott_nonoro)
286	write(,)" calllott_nonoro = ",calllott_nonoro
287
288	! rocket dust storm injection scheme
289	write(,)"call rocket dust storm and slope lifting",
290	& " parametrization"
291	rdstorm=.false. ! default value
292	call getin("rdstorm",rdstorm)
293	write(,)" rdstorm = ",rdstorm
294	! rocket dust storm detrainment coefficient
295	coeff_detrainment=0. ! default value
296	call getin("coeff_detrainment",coeff_detrainment)
297	write(,)" coeff_detrainment = ",coeff_detrainment
298
299	write(,)"rad.transfer is computed every iradia",
300	& " physical timestep"
301	iradia=1 ! default value
302	call getin("iradia",iradia)
303	write(,)" iradia = ",iradia
304
305
306	write(,)"Output of the exchange coefficient mattrix ?",
307	& "(for diagnostics only)"
308	callg2d=.false. ! default value
309	call getin("callg2d",callg2d)
310	write(,)" callg2d = ",callg2d
311
312	write(,)"Rayleigh scattering : (should be .false. for now)"
313	rayleigh=.false.
314	call getin("rayleigh",rayleigh)
315	write(,)" rayleigh = ",rayleigh
316
317
318	! TRACERS:
319
320	! dustbin
321	write(,)"Transported dust ? (if >0, use 'dustbin' dust bins)"
322	dustbin=0 ! default value
323	call getin("dustbin",dustbin)
324	write(,)" dustbin = ",dustbin
325	! active
326	write(,)"Radiatively active dust ? (matters if dustbin>0)"
327	active=.false. ! default value
328	call getin("active",active)
329	write(,)" active = ",active
330
331	! Test of incompatibility:
332	! if active is used, then dustbin should be > 0
333
334	if (active.and.(dustbin.lt.1)) then
335	print*,'if active is used, then dustbin should > 0'
336	stop
337	endif
338	! doubleq
339	write(,)"use mass and number mixing ratios to predict",
340	& " dust size ?"
341	doubleq=.false. ! default value
342	call getin("doubleq",doubleq)
343	write(,)" doubleq = ",doubleq
344	! submicron
345	submicron=.false. ! default value
346	call getin("submicron",submicron)
347	write(,)" submicron = ",submicron
348
349	! Test of incompatibility:
350	! if doubleq is used, then dustbin should be 2
351
352	if (doubleq.and.(dustbin.ne.2)) then
353	print*,'if doubleq is used, then dustbin should be 2'
354	stop
355	endif
356	if (doubleq.and.submicron.and.(nq.LT.3)) then
357	print*,'If doubleq is used with a submicron tracer,'
358	print*,' then the number of tracers has to be'
359	print*,' larger than 3.'
360	stop
361	endif
362
363	! lifting
364	write(,)"dust lifted by GCM surface winds ?"
365	lifting=.false. ! default value
366	call getin("lifting",lifting)
367	write(,)" lifting = ",lifting
368
369	! Test of incompatibility:
370	! if lifting is used, then dustbin should be > 0
371
372	if (lifting.and.(dustbin.lt.1)) then
373	print*,'if lifting is used, then dustbin should > 0'
374	stop
375	endif
376
377	! dust injection scheme
378	dustinjection=0 ! default: no injection scheme
379	call getin("dustinjection",dustinjection)
380	write(,)" dustinjection = ",dustinjection
381	! dust injection scheme coefficient
382	coeff_injection=1. ! default value
383	call getin("coeff_injection",coeff_injection)
384	write(,)" coeff_in,jection = ",coeff_injection
385	! timing for dust injection
386	ti_injection=10. ! default value
387	tf_injection=12. ! default value
388	call getin("ti_injection",ti_injection)
389	write(,)" ti_injection = ",ti_injection
390	call getin("tf_injection",tf_injection)
391	write(,)" tf_injection = ",tf_injection
392
393	! free evolving dust
394	! freedust=true just says that there is no lifting and no dust opacity scaling.
395	write(,)"dust lifted by GCM surface winds ?"
396	freedust=.false. ! default value
397	call getin("freedust",freedust)
398	write(,)" freedust = ",freedust
399	if (freedust.and..not.doubleq) then
400	print*,'freedust should be used with doubleq !'
401	stop
402	endif
403	#ifndef MESOSCALE
404	! this test is valid in GCM case
405	! ... not in mesoscale case, for we want to activate mesoscale lifting
406	if (freedust.and.dustinjection.eq.0)then
407	if(lifting) then
408	print*,'if freedust is used and dustinjection = 0,
409	& then lifting should not be used'
410	stop
411	endif
412	endif
413	#endif
414	if (dustinjection.eq.1)then
415	if(.not.lifting) then
416	print*,"if dustinjection=1, then lifting should be true"
417	stop
418	endif
419	if(.not.freedust) then
420	print*,"if dustinjection=1, then freedust should be true"
421	stop
422	endif
423	endif
424	! rocket dust storm
425	! Test of incompatibility:
426	! if rdstorm is used, then doubleq should be true
427	if (rdstorm.and..not.doubleq) then
428	print*,'if rdstorm is used, then doubleq should be used !'
429	stop
430	endif
431	if (rdstorm.and..not.active) then
432	print*,'if rdstorm is used, then active should be used !'
433	stop
434	endif
435	if (rdstorm.and..not.lifting) then
436	print*,'if rdstorm is used, then lifting should be used !'
437	stop
438	endif
439	if (rdstorm.and..not.freedust) then
440	print*,'if rdstorm is used, then freedust should be used !'
441	stop
442	endif
443	if (rdstorm.and.(dustinjection.eq.0)) then
444	print*,'if rdstorm is used, then dustinjection should
445	& be used !'
446	stop
447	endif
448	! Dust IR opacity
449	write(,)" Wavelength for infrared opacity of dust ?"
450	write(,)" Choices are:"
451	write(,)" tes --- > 9.3 microns [default]"
452	write(,)" mcs --- > 21.6 microns"
453	!
454	! WARNING WARNING WARNING WARNING WARNING WARNING
455	!
456	! BEFORE ADDING A NEW VALUE, BE SURE THAT THE
457	! CORRESPONDING WAVELENGTH IS IN THE LOOKUP TABLE,
458	! OR AT LEAST NO TO FAR, TO AVOID FALLACIOUS INTERPOLATIONS.
459	!
460	dustiropacity="tes" !default value - is expected to shift to mcs one day
461	call getin("dustiropacity",dustiropacity)
462	write(,)" dustiropacity = ",trim(dustiropacity)
463	select case (trim(dustiropacity))
464	case ("tes")
465	dustrefir = 9.3E-6
466	case ("mcs")
467	dustrefir = 21.6E-6
468	case default
469	write(,) trim(dustiropacity),
470	& " is not a valid option for dustiropacity"
471	stop
472	end select
473
474	! callddevil
475	write(,)" dust lifted by dust devils ?"
476	callddevil=.false. !default value
477	call getin("callddevil",callddevil)
478	write(,)" callddevil = ",callddevil
479
480	! Test of incompatibility:
481	! if dustdevil is used, then dustbin should be > 0
482
483	if (callddevil.and.(dustbin.lt.1)) then
484	print*,'if dustdevil is used, then dustbin should > 0'
485	stop
486	endif
487	! sedimentation
488	write(,) "Gravitationnal sedimentation ?"
489	sedimentation=.true. ! default value
490	call getin("sedimentation",sedimentation)
491	write(,) " sedimentation = ",sedimentation
492	! activice
493	write(,) "Radiatively active transported atmospheric ",
494	& "water ice ?"
495	activice=.false. ! default value
496	call getin("activice",activice)
497	write(,) " activice = ",activice
498	! water
499	write(,) "Compute water cycle ?"
500	water=.false. ! default value
501	call getin("water",water)
502	write(,) " water = ",water
503
504	! sub-grid cloud fraction: fixed
505	write(,) "Fixed cloud fraction?"
506	CLFfixval=1.0 ! default value
507	call getin("CLFfixval",CLFfixval)
508	write(,) "CLFfixval=",CLFfixval
509	! sub-grid cloud fraction: varying
510	write(,) "Use partial nebulosity?"
511	CLFvarying=.false. ! default value
512	call getin("CLFvarying",CLFvarying)
513	write(,)"CLFvarying=",CLFvarying
514
515	!CO2 clouds scheme?
516	write(,) "Compute CO2 clouds (implies microphysical scheme)?"
517	co2clouds=.false. ! default value
518	call getin("co2clouds",co2clouds)
519	write(,) " co2clouds = ",co2clouds
520	!Can water ice particles serve as CCN for CO2clouds
521	write(,) "Use water ice as CO2 clouds CCN ?"
522	co2useh2o=.false. ! default value
523	call getin("co2useh2o",co2useh2o)
524	write(,) " co2useh2o = ",co2useh2o
525	!Do we allow a supply of meteoritic paricles to serve as CO2 ice CCN?
526	write(,) "Supply meteoritic particle for CO2 clouds ?"
527	meteo_flux=.false. !Default value
528	call getin("meteo_flux",meteo_flux)
529	write(,) " meteo_flux = ",meteo_flux
530	!Do we allow a sub-grid temperature distribution for the CO2 microphysics
531	write(,) "sub-grid temperature distribution for CO2 clouds?"
532	CLFvaryingCO2=.false. !Default value
533	call getin("CLFvaryingCO2",CLFvaryingCO2)
534	write(,) " CLFvaryingCO2 = ",CLFvaryingCO2
535	!Amplitude of the sub-grid temperature distribution for the CO2 microphysics
536	write(,) "sub-grid temperature amplitude for CO2 clouds?"
537	spantCO2=0 !Default value
538	call getin("spantCO2",spantCO2)
539	write(,) " spantCO2 = ",spantCO2
540	!Do you want to filter the sub-grid T distribution by a Saturation index?
541	write(,) "filter sub-grid temperature by Saturation index?"
542	satindexco2=.true.
543	call getin("satindexco2",satindexco2)
544	write(,) " satindexco2 = ",satindexco2
545
546
547	! thermal inertia feedback
548	write(,) "Activate the thermal inertia feedback ?"
549	tifeedback=.false. ! default value
550	call getin("tifeedback",tifeedback)
551	write(,) " tifeedback = ",tifeedback
552
553	! Test of incompatibility:
554
555	if (tifeedback.and..not.water) then
556	print*,'if tifeedback is used,'
557	print*,'water should be used too'
558	stop
559	endif
560
561	if (tifeedback.and..not.callsoil) then
562	print*,'if tifeedback is used,'
563	print*,'callsoil should be used too'
564	stop
565	endif
566
567	if (activice.and..not.water) then
568	print*,'if activice is used, water should be used too'
569	stop
570	endif
571
572	if (water.and..not.tracer) then
573	print*,'if water is used, tracer should be used too'
574	stop
575	endif
576
577	! water ice clouds effective variance distribution for sedimentaion
578	write(,) "Sed effective variance for water ice clouds ?"
579	nuice_sed=0.45
580	call getin("nuice_sed",nuice_sed)
581	write(,) "water_param nueff Sedimentation:", nuice_sed
582
583	write(,) "Sed effective variance for CO2 clouds ?"
584	nuiceco2_sed=0.45
585	call getin("nuiceco2_sed",nuiceco2_sed)
586	write(,) "CO2 nueff Sedimentation:", nuiceco2_sed
587
588	write(,) "REF effective variance for CO2 clouds ?"
589	nuiceco2_ref=0.45
590	call getin("nuiceco2_ref",nuiceco2_ref)
591	write(,) "CO2 nueff Sedimentation:", nuiceco2_ref
592
593	write(,) "REF effective variance for water clouds ?"
594	nuice_ref=0.45
595	call getin("nuice_ref",nuice_ref)
596	write(,) "CO2 nueff Sedimentation:", nuice_ref
597
598
599	! ccn factor if no scavenging
600	write(,) "water param CCN reduc. factor ?"
601	ccn_factor = 4.5
602	call getin("ccn_factor",ccn_factor)
603	write(,)" ccn_factor = ",ccn_factor
604	write(,)"Careful: only used when microphys=F, otherwise"
605	write(,)"the contact parameter is used instead;"
606
607	! microphys
608	write(,)"Microphysical scheme for water-ice clouds?"
609	microphys=.false. ! default value
610	call getin("microphys",microphys)
611	write(,)" microphys = ",microphys
612
613	! supersat
614	write(,)"Allow super-saturation of water vapor?"
615	supersat=.true. ! default value
616	call getin("supersat",supersat)
617	write(,)"supersat = ",supersat
618
619	! microphysical parameter contact
620	write(,) "water contact parameter ?"
621	mteta = 0.95
622	call getin("mteta",mteta)
623	write(,) "mteta = ", mteta
624
625	! scavenging
626	write(,)"Dust scavenging by H2O/CO2 snowfall ?"
627	scavenging=.false. ! default value
628	call getin("scavenging",scavenging)
629	write(,)" scavenging = ",scavenging
630
631
632	! Test of incompatibility:
633	! if scavenging is used, then dustbin should be > 0
634
635	if ((microphys.and..not.doubleq).or.
636	& (microphys.and..not.water)) then
637	print*,'if microphys is used, then doubleq,'
638	print*,'and water must be used!'
639	stop
640	endif
641	if (microphys.and..not.scavenging) then
642	print*,''
643	print*,'----------------WARNING-----------------'
644	print*,'microphys is used without scavenging !!!'
645	print*,'----------------WARNING-----------------'
646	print*,''
647	endif
648
649	if ((scavenging.and..not.microphys).or.
650	& (scavenging.and.(dustbin.lt.1)))then
651	print*,'if scavenging is used, then microphys'
652	print*,'must be used!'
653	stop
654	endif
655
656	! Test of incompatibility:
657
658	write(,) "Permanent water caps at poles ?",
659	& " .true. is RECOMMENDED"
660	write(,) "(with .true., North cap is a source of water ",
661	& "and South pole is a cold trap)"
662	caps=.true. ! default value
663	call getin("caps",caps)
664	write(,) " caps = ",caps
665
666	! albedo_h2o_ice
667	write(,) "water ice albedo ?"
668	albedo_h2o_ice=0.45
669	call getin("albedo_h2o_ice",albedo_h2o_ice)
670	write(,) " albedo_h2o_ice = ",albedo_h2o_ice
671	! inert_h2o_ice
672	write(,) "water ice thermal inertia ?"
673	inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
674	call getin("inert_h2o_ice",inert_h2o_ice)
675	write(,) " inert_h2o_ice = ",inert_h2o_ice
676	! frost_albedo_threshold
677	write(,) "frost thickness threshold for albedo ?"
678	frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2)
679	call getin("frost_albedo_threshold",
680	& frost_albedo_threshold)
681	write(,) " frost_albedo_threshold = ",
682	& frost_albedo_threshold
683
684	! call Titus crocus line -- DEFAULT IS NONE
685	write(,) "Titus crocus line ?"
686	tituscap=.false. ! default value
687	call getin("tituscap",tituscap)
688	write(,) "tituscap",tituscap
689
690
691	write(,) "photochemistry: include chemical species"
692	photochem=.false. ! default value
693	call getin("photochem",photochem)
694	write(,) " photochem = ",photochem
695
696
697	! SCATTERERS
698	write(,) "how many scatterers?"
699	naerkind=1 ! default value
700	call getin("naerkind",naerkind)
701	write(,)" naerkind = ",naerkind
702
703	! Test of incompatibility
704	c Logical tests for radiatively active water-ice clouds:
705	IF ( (activice.AND.(.NOT.water)).OR.
706	& (activice.AND.(naerkind.LT.2)) ) THEN
707	WRITE(,) 'If activice is TRUE, water has to be set'
708	WRITE(,) 'to TRUE, and "naerkind" must be at least'
709	WRITE(,) 'equal to 2.'
710	CALL ABORT
711	ENDIF
712
713	!------------------------------------------
714	!------------------------------------------
715	! once naerkind is known allocate arrays
716	! -- we do it here and not in phys_var_init
717	! -- because we need to know naerkind first
718	CALL ini_scatterers(ngrid,nlayer)
719	!------------------------------------------
720	!------------------------------------------
721
722
723	c Please name the different scatterers here ----------------
724	name_iaer(1) = "dust_conrath" !! default choice is good old Conrath profile
725	IF (doubleq.AND.active) name_iaer(1) = "dust_doubleq" !! two-moment scheme
726	if (nq.gt.1) then
727	! trick to avoid problems compiling with 1 tracer
728	! and picky compilers who know name_iaer(2) is out of bounds
729	j=2
730	IF (rdstorm.AND..NOT.activice) name_iaer(2) =
731	& "stormdust_doubleq" !! storm dust two-moment scheme
732	IF (rdstorm.AND.water.AND.activice) name_iaer(3) =
733	& "stormdust_doubleq"
734	IF (water.AND.activice) name_iaer(j) = "h2o_ice" !! radiatively-active clouds
735	IF (submicron.AND.active) name_iaer(j) = "dust_submicron" !! JBM experimental stuff
736	endif ! of if (nq.gt.1)
737
738	c ----------------------------------------------------------
739
740	! THERMOSPHERE
741
742	write(,) "call thermosphere ?"
743	callthermos=.false. ! default value
744	call getin("callthermos",callthermos)
745	write(,) " callthermos = ",callthermos
746
747
748	write(,) " water included without cycle ",
749	& "(only if water=.false.)"
750	thermoswater=.false. ! default value
751	call getin("thermoswater",thermoswater)
752	write(,) " thermoswater = ",thermoswater
753
754	write(,) "call thermal conduction ?",
755	& " (only if callthermos=.true.)"
756	callconduct=.false. ! default value
757	call getin("callconduct",callconduct)
758	write(,) " callconduct = ",callconduct
759
760	write(,) "call EUV heating ?",
761	& " (only if callthermos=.true.)"
762	calleuv=.false. ! default value
763	call getin("calleuv",calleuv)
764	write(,) " calleuv = ",calleuv
765
766	write(,) "call molecular viscosity ?",
767	& " (only if callthermos=.true.)"
768	callmolvis=.false. ! default value
769	call getin("callmolvis",callmolvis)
770	write(,) " callmolvis = ",callmolvis
771
772	write(,) "call molecular diffusion ?",
773	& " (only if callthermos=.true.)"
774	callmoldiff=.false. ! default value
775	call getin("callmoldiff",callmoldiff)
776	write(,) " callmoldiff = ",callmoldiff
777
778
779	write(,) "call thermospheric photochemistry ?",
780	& " (only if callthermos=.true.)"
781	thermochem=.false. ! default value
782	call getin("thermochem",thermochem)
783	write(,) " thermochem = ",thermochem
784
785	write(,) "Method to include solar variability"
786	write(,) "0-> fixed value of E10.7 (fixed_euv_value); ",
787	& "1-> daily evolution of E10.7 (for given solvaryear)"
788	solvarmod=1
789	call getin("solvarmod",solvarmod)
790	write(,) " solvarmod = ",solvarmod
791
792	write(,) "Fixed euv (for solvarmod==0) 10.7 value?"
793	write(,) " (min=80 , ave=140, max=320)"
794	fixed_euv_value=140 ! default value
795	call getin("fixed_euv_value",fixed_euv_value)
796	write(,) " fixed_euv_value = ",fixed_euv_value
797
798	write(,) "Solar variability as observed for MY: "
799	write(,) "Only if solvarmod=1"
800	solvaryear=24
801	call getin("solvaryear",solvaryear)
802	write(,) " solvaryear = ",solvaryear
803
804	write(,) "UV heating efficiency:",
805	& "measured values between 0.19 and 0.23 (Fox et al. 1996)",
806	& "lower values may be used to compensate low 15 um cooling"
807	euveff=0.21 !default value
808	call getin("euveff",euveff)
809	write(,) " euveff = ", euveff
810
811
812	if (.not.callthermos) then
813	if (thermoswater) then
814	print*,'if thermoswater is set, callthermos must be true'
815	stop
816	endif
817	if (callconduct) then
818	print*,'if callconduct is set, callthermos must be true'
819	stop
820	endif
821	if (calleuv) then
822	print*,'if calleuv is set, callthermos must be true'
823	stop
824	endif
825	if (callmolvis) then
826	print*,'if callmolvis is set, callthermos must be true'
827	stop
828	endif
829	if (callmoldiff) then
830	print*,'if callmoldiff is set, callthermos must be true'
831	stop
832	endif
833	if (thermochem) then
834	print*,'if thermochem is set, callthermos must be true'
835	stop
836	endif
837	endif
838
839	! Test of incompatibility:
840	! if photochem is used, then water should be used too
841
842	if (photochem.and..not.water) then
843	print*,'if photochem is used, water should be used too'
844	stop
845	endif
846
847	! if callthermos is used, then thermoswater should be used too
848	! (if water not used already)
849
850	if (callthermos .and. .not.water) then
851	if (callthermos .and. .not.thermoswater) then
852	print*,'if callthermos is used, water or thermoswater
853	& should be used too'
854	stop
855	endif
856	endif
857
858	PRINT*,'--------------------------------------------'
859	PRINT*
860	PRINT*
861	ELSE
862	write(,)
863	write(,) 'Cannot read file callphys.def. Is it here ?'
864	stop
865	ENDIF
866
867	8000 FORMAT(t5,a12,l8)
868	8001 FORMAT(t5,a12,i8)
869
870	PRINT*
871	PRINT*,'conf_phys: daysec',daysec
872	PRINT*
873	PRINT*,'conf_phys: The radiative transfer is computed:'
874	PRINT*,' each ',iradia,' physical time-step'
875	PRINT,' or each ',iradiadtphys,' seconds'
876	PRINT*
877	! --------------------------------------------------------------
878	! Managing the Longwave radiative transfer
879	! --------------------------------------------------------------
880
881	! In most cases, the run just use the following values :
882	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
883	callemis=.true.
884	! ilwd=10*int(daysec/dtphys) ! bug before 22/10/01
885	ilwd=1
886	ilwn=1 !2
887	ilwb=1 !2
888	linear=.true.
889	ncouche=3
890	alphan=0.4
891	semi=0
892
893	! BUT people working hard on the LW may want to read them in 'radia.def'
894	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
895
896	OPEN(99,file='radia.def',status='old',form='formatted'
897	. ,iostat=ierr)
898	IF(ierr.EQ.0) THEN
899	write(,) 'conf_phys: Reading radia.def !!!'
900	READ(99,fmt='(a)') ch1
901	READ(99,*) callemis
902	WRITE(*,8000) ch1,callemis
903
904	READ(99,fmt='(a)') ch1
905	READ(99,*) iradia
906	WRITE(*,8001) ch1,iradia
907
908	READ(99,fmt='(a)') ch1
909	READ(99,*) ilwd
910	WRITE(*,8001) ch1,ilwd
911
912	READ(99,fmt='(a)') ch1
913	READ(99,*) ilwn
914	WRITE(*,8001) ch1,ilwn
915
916	READ(99,fmt='(a)') ch1
917	READ(99,*) linear
918	WRITE(*,8000) ch1,linear
919
920	READ(99,fmt='(a)') ch1
921	READ(99,*) ncouche
922	WRITE(*,8001) ch1,ncouche
923
924	READ(99,fmt='(a)') ch1
925	READ(99,*) alphan
926	WRITE(,) ch1,alphan
927
928	READ(99,fmt='(a)') ch1
929	READ(99,*) ilwb
930	WRITE(*,8001) ch1,ilwb
931
932
933	READ(99,fmt='(a)') ch1
934	READ(99,'(l1)') callg2d
935	WRITE(*,8000) ch1,callg2d
936
937	READ(99,fmt='(a)') ch1
938	READ(99,*) semi
939	WRITE(,) ch1,semi
940	end if
941	CLOSE(99)
942
943	END

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

Download in other formats:

Original Format