Context Navigation

conf_phys.F @ 2184

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

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

Download in other formats:

Original Format