Context Navigation

conf_phys.F @ 2266

Last change on this file since 2266 was 2218, checked in by emillour, 5 years ago
Mars GCM: Change "latentheat" flag to a more descriptive "latentheat_surfwater" and set its default value to .true. EM
File size: 36.1 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	! entrainment by slope wind scheme
301	write(,)"call slope wind lifting parametrization"
302	slpwind=.false. ! default value
303	call getin("slpwind",slpwind)
304	write(,)" slpwind = ",slpwind
305
306	! latent heat release from ground water ice sublimation/condensation
307	write(,)"latent heat release during sublimation",
308	& " /condensation of ground water ice"
309	latentheat_surfwater=.true. ! default value
310	call getin("latentheat_surfwater",latentheat_surfwater)
311	write(,)" latentheat_surfwater = ",latentheat_surfwater
312
313	write(,)"rad.transfer is computed every iradia",
314	& " physical timestep"
315	iradia=1 ! default value
316	call getin("iradia",iradia)
317	write(,)" iradia = ",iradia
318
319
320	write(,)"Output of the exchange coefficient mattrix ?",
321	& "(for diagnostics only)"
322	callg2d=.false. ! default value
323	call getin("callg2d",callg2d)
324	write(,)" callg2d = ",callg2d
325
326	write(,)"Rayleigh scattering : (should be .false. for now)"
327	rayleigh=.false.
328	call getin("rayleigh",rayleigh)
329	write(,)" rayleigh = ",rayleigh
330
331
332	! TRACERS:
333
334	! dustbin
335	write(,)"Transported dust ? (if >0, use 'dustbin' dust bins)"
336	dustbin=0 ! default value
337	call getin("dustbin",dustbin)
338	write(,)" dustbin = ",dustbin
339	! active
340	write(,)"Radiatively active dust ? (matters if dustbin>0)"
341	active=.false. ! default value
342	call getin("active",active)
343	write(,)" active = ",active
344
345	! Test of incompatibility:
346	! if active is used, then dustbin should be > 0
347
348	if (active.and.(dustbin.lt.1)) then
349	print*,'if active is used, then dustbin should > 0'
350	stop
351	endif
352	! doubleq
353	write(,)"use mass and number mixing ratios to predict",
354	& " dust size ?"
355	doubleq=.false. ! default value
356	call getin("doubleq",doubleq)
357	write(,)" doubleq = ",doubleq
358	! submicron
359	submicron=.false. ! default value
360	call getin("submicron",submicron)
361	write(,)" submicron = ",submicron
362
363	! Test of incompatibility:
364	! if doubleq is used, then dustbin should be 2
365
366	if (doubleq.and.(dustbin.ne.2)) then
367	print*,'if doubleq is used, then dustbin should be 2'
368	stop
369	endif
370	if (doubleq.and.submicron.and.(nq.LT.3)) then
371	print*,'If doubleq is used with a submicron tracer,'
372	print*,' then the number of tracers has to be'
373	print*,' larger than 3.'
374	stop
375	endif
376
377	! lifting
378	write(,)"dust lifted by GCM surface winds ?"
379	lifting=.false. ! default value
380	call getin("lifting",lifting)
381	write(,)" lifting = ",lifting
382
383	! Test of incompatibility:
384	! if lifting is used, then dustbin should be > 0
385
386	if (lifting.and.(dustbin.lt.1)) then
387	print*,'if lifting is used, then dustbin should > 0'
388	stop
389	endif
390
391	! dust injection scheme
392	dustinjection=0 ! default: no injection scheme
393	call getin("dustinjection",dustinjection)
394	write(,)" dustinjection = ",dustinjection
395	! dust injection scheme coefficient
396	coeff_injection=1. ! default value
397	call getin("coeff_injection",coeff_injection)
398	write(,)" coeff_in,jection = ",coeff_injection
399	! timing for dust injection
400	ti_injection=10. ! default value
401	tf_injection=12. ! default value
402	call getin("ti_injection",ti_injection)
403	write(,)" ti_injection = ",ti_injection
404	call getin("tf_injection",tf_injection)
405	write(,)" tf_injection = ",tf_injection
406
407	! free evolving dust
408	! freedust=true just says that there is no lifting and no dust opacity scaling.
409	write(,)"dust lifted by GCM surface winds ?"
410	freedust=.false. ! default value
411	call getin("freedust",freedust)
412	write(,)" freedust = ",freedust
413	if (freedust.and..not.doubleq) then
414	print*,'freedust should be used with doubleq !'
415	stop
416	endif
417	#ifndef MESOSCALE
418	! this test is valid in GCM case
419	! ... not in mesoscale case, for we want to activate mesoscale lifting
420	if (freedust.and.dustinjection.eq.0)then
421	if(lifting) then
422	print*,'if freedust is used and dustinjection = 0,
423	& then lifting should not be used'
424	stop
425	endif
426	endif
427	#endif
428	if (dustinjection.eq.1)then
429	if(.not.lifting) then
430	print*,"if dustinjection=1, then lifting should be true"
431	stop
432	endif
433	if(.not.freedust) then
434	print*,"if dustinjection=1, then freedust should be true"
435	stop
436	endif
437	endif
438	! rocket dust storm and entrainment by slope wind
439	! Test of incompatibility:
440	! if rdstorm or slpwind is used, then doubleq should be true
441	if ((rdstorm.or.slpwind).and..not.doubleq) then
442	print*,'if rdstorm or slpwind is used, then doubleq
443	& should be used !'
444	stop
445	endif
446	if ((rdstorm.or.slpwind).and..not.active) then
447	print*,'if rdstorm or slpwind is used, then active
448	& should be used !'
449	stop
450	endif
451	if (rdstorm.and..not.lifting) then
452	print*,'if rdstorm is used, then lifting
453	& should be used !'
454	stop
455	endif
456	if ((rdstorm.or.slpwind).and..not.freedust) then
457	print*,'if rdstorm or slpwind is used, then freedust
458	& should be used !'
459	stop
460	endif
461	if (rdstorm.and.(dustinjection.eq.0)) then
462	print*,'if rdstorm is used, then dustinjection
463	& should be used !'
464	stop
465	endif
466	! Dust IR opacity
467	write(,)" Wavelength for infrared opacity of dust ?"
468	write(,)" Choices are:"
469	write(,)" tes --- > 9.3 microns [default]"
470	write(,)" mcs --- > 21.6 microns"
471	!
472	! WARNING WARNING WARNING WARNING WARNING WARNING
473	!
474	! BEFORE ADDING A NEW VALUE, BE SURE THAT THE
475	! CORRESPONDING WAVELENGTH IS IN THE LOOKUP TABLE,
476	! OR AT LEAST NO TO FAR, TO AVOID FALLACIOUS INTERPOLATIONS.
477	!
478	dustiropacity="tes" !default value - is expected to shift to mcs one day
479	call getin("dustiropacity",dustiropacity)
480	write(,)" dustiropacity = ",trim(dustiropacity)
481	select case (trim(dustiropacity))
482	case ("tes")
483	dustrefir = 9.3E-6
484	case ("mcs")
485	dustrefir = 21.6E-6
486	case default
487	write(,) trim(dustiropacity),
488	& " is not a valid option for dustiropacity"
489	stop
490	end select
491
492	! callddevil
493	write(,)" dust lifted by dust devils ?"
494	callddevil=.false. !default value
495	call getin("callddevil",callddevil)
496	write(,)" callddevil = ",callddevil
497
498	! Test of incompatibility:
499	! if dustdevil is used, then dustbin should be > 0
500
501	if (callddevil.and.(dustbin.lt.1)) then
502	print*,'if dustdevil is used, then dustbin should > 0'
503	stop
504	endif
505	! sedimentation
506	write(,) "Gravitationnal sedimentation ?"
507	sedimentation=.true. ! default value
508	call getin("sedimentation",sedimentation)
509	write(,) " sedimentation = ",sedimentation
510	! activice
511	write(,) "Radiatively active transported atmospheric ",
512	& "water ice ?"
513	activice=.false. ! default value
514	call getin("activice",activice)
515	write(,) " activice = ",activice
516	! water
517	write(,) "Compute water cycle ?"
518	water=.false. ! default value
519	call getin("water",water)
520	write(,) " water = ",water
521
522	! sub-grid cloud fraction: fixed
523	write(,) "Fixed cloud fraction?"
524	CLFfixval=1.0 ! default value
525	call getin("CLFfixval",CLFfixval)
526	write(,) "CLFfixval=",CLFfixval
527	! sub-grid cloud fraction: varying
528	write(,) "Use partial nebulosity?"
529	CLFvarying=.false. ! default value
530	call getin("CLFvarying",CLFvarying)
531	write(,)"CLFvarying=",CLFvarying
532
533	!CO2 clouds scheme?
534	write(,) "Compute CO2 clouds (implies microphysical scheme)?"
535	co2clouds=.false. ! default value
536	call getin("co2clouds",co2clouds)
537	write(,) " co2clouds = ",co2clouds
538	!Can water ice particles serve as CCN for CO2clouds
539	write(,) "Use water ice as CO2 clouds CCN ?"
540	co2useh2o=.false. ! default value
541	call getin("co2useh2o",co2useh2o)
542	write(,) " co2useh2o = ",co2useh2o
543	!Do we allow a supply of meteoritic paricles to serve as CO2 ice CCN?
544	write(,) "Supply meteoritic particle for CO2 clouds ?"
545	meteo_flux=.false. !Default value
546	call getin("meteo_flux",meteo_flux)
547	write(,) " meteo_flux = ",meteo_flux
548	!Do we allow a sub-grid temperature distribution for the CO2 microphysics
549	write(,) "sub-grid temperature distribution for CO2 clouds?"
550	CLFvaryingCO2=.false. !Default value
551	call getin("CLFvaryingCO2",CLFvaryingCO2)
552	write(,) " CLFvaryingCO2 = ",CLFvaryingCO2
553	!Amplitude of the sub-grid temperature distribution for the CO2 microphysics
554	write(,) "sub-grid temperature amplitude for CO2 clouds?"
555	spantCO2=0 !Default value
556	call getin("spantCO2",spantCO2)
557	write(,) " spantCO2 = ",spantCO2
558	!Do you want to filter the sub-grid T distribution by a Saturation index?
559	write(,) "filter sub-grid temperature by Saturation index?"
560	satindexco2=.true.
561	call getin("satindexco2",satindexco2)
562	write(,) " satindexco2 = ",satindexco2
563
564
565	! thermal inertia feedback
566	write(,) "Activate the thermal inertia feedback ?"
567	tifeedback=.false. ! default value
568	call getin("tifeedback",tifeedback)
569	write(,) " tifeedback = ",tifeedback
570
571	! Test of incompatibility:
572
573	if (tifeedback.and..not.water) then
574	print*,'if tifeedback is used,'
575	print*,'water should be used too'
576	stop
577	endif
578
579	if (tifeedback.and..not.callsoil) then
580	print*,'if tifeedback is used,'
581	print*,'callsoil should be used too'
582	stop
583	endif
584
585	if (activice.and..not.water) then
586	print*,'if activice is used, water should be used too'
587	stop
588	endif
589
590	if (water.and..not.tracer) then
591	print*,'if water is used, tracer should be used too'
592	stop
593	endif
594
595	! water ice clouds effective variance distribution for sedimentaion
596	write(,) "Sed effective variance for water ice clouds ?"
597	nuice_sed=0.45
598	call getin("nuice_sed",nuice_sed)
599	write(,) "water_param nueff Sedimentation:", nuice_sed
600
601	write(,) "Sed effective variance for CO2 clouds ?"
602	nuiceco2_sed=0.45
603	call getin("nuiceco2_sed",nuiceco2_sed)
604	write(,) "CO2 nueff Sedimentation:", nuiceco2_sed
605
606	write(,) "REF effective variance for CO2 clouds ?"
607	nuiceco2_ref=0.45
608	call getin("nuiceco2_ref",nuiceco2_ref)
609	write(,) "CO2 nueff Sedimentation:", nuiceco2_ref
610
611	write(,) "REF effective variance for water clouds ?"
612	nuice_ref=0.45
613	call getin("nuice_ref",nuice_ref)
614	write(,) "CO2 nueff Sedimentation:", nuice_ref
615
616
617	! ccn factor if no scavenging
618	write(,) "water param CCN reduc. factor ?"
619	ccn_factor = 4.5
620	call getin("ccn_factor",ccn_factor)
621	write(,)" ccn_factor = ",ccn_factor
622	write(,)"Careful: only used when microphys=F, otherwise"
623	write(,)"the contact parameter is used instead;"
624
625	! microphys
626	write(,)"Microphysical scheme for water-ice clouds?"
627	microphys=.false. ! default value
628	call getin("microphys",microphys)
629	write(,)" microphys = ",microphys
630
631	! supersat
632	write(,)"Allow super-saturation of water vapor?"
633	supersat=.true. ! default value
634	call getin("supersat",supersat)
635	write(,)"supersat = ",supersat
636
637	! microphysical parameter contact
638	write(,) "water contact parameter ?"
639	mteta = 0.95
640	call getin("mteta",mteta)
641	write(,) "mteta = ", mteta
642
643	! scavenging
644	write(,)"Dust scavenging by H2O/CO2 snowfall ?"
645	scavenging=.false. ! default value
646	call getin("scavenging",scavenging)
647	write(,)" scavenging = ",scavenging
648
649
650	! Test of incompatibility:
651	! if scavenging is used, then dustbin should be > 0
652
653	if ((microphys.and..not.doubleq).or.
654	& (microphys.and..not.water)) then
655	print*,'if microphys is used, then doubleq,'
656	print*,'and water must be used!'
657	stop
658	endif
659	if (microphys.and..not.scavenging) then
660	print*,''
661	print*,'----------------WARNING-----------------'
662	print*,'microphys is used without scavenging !!!'
663	print*,'----------------WARNING-----------------'
664	print*,''
665	endif
666
667	if ((scavenging.and..not.microphys).or.
668	& (scavenging.and.(dustbin.lt.1)))then
669	print*,'if scavenging is used, then microphys'
670	print*,'must be used!'
671	stop
672	endif
673
674	! Instantaneous scavenging by CO2
675	! -> expected to be replaced by scavenging with microphysics (flag scavenging) one day
676	write(,)"Dust scavenging by instantaneous CO2 snowfall ?"
677	scavco2cond=.false. ! default value
678	call getin("scavco2cond",scavco2cond)
679	write(,)" scavco2cond = ",scavco2cond
680	! Test of incompatibility:
681	! if scavco2cond is used, then dustbin should be > 0
682	if (scavco2cond.and.(dustbin.lt.1))then
683	print*,'if scavco2cond is used, then dustbin should be > 0'
684	stop
685	endif
686	! if co2clouds is used, then there is no need for scavco2cond
687	if (co2clouds.and.scavco2cond) then
688	print*,''
689	print*,'----------------WARNING-----------------'
690	print*,' microphys scavenging is used so '
691	print*,' no need for scavco2cond !!! '
692	print*,'----------------WARNING-----------------'
693	print*,''
694	stop
695	endif
696
697	! Test of incompatibility:
698
699	write(,) "Permanent water caps at poles ?",
700	& " .true. is RECOMMENDED"
701	write(,) "(with .true., North cap is a source of water ",
702	& "and South pole is a cold trap)"
703	caps=.true. ! default value
704	call getin("caps",caps)
705	write(,) " caps = ",caps
706
707	! albedo_h2o_ice
708	write(,) "water ice albedo ?"
709	albedo_h2o_ice=0.45
710	call getin("albedo_h2o_ice",albedo_h2o_ice)
711	write(,) " albedo_h2o_ice = ",albedo_h2o_ice
712	! inert_h2o_ice
713	write(,) "water ice thermal inertia ?"
714	inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
715	call getin("inert_h2o_ice",inert_h2o_ice)
716	write(,) " inert_h2o_ice = ",inert_h2o_ice
717	! frost_albedo_threshold
718	write(,) "frost thickness threshold for albedo ?"
719	frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2)
720	call getin("frost_albedo_threshold",
721	& frost_albedo_threshold)
722	write(,) " frost_albedo_threshold = ",
723	& frost_albedo_threshold
724
725	! call Titus crocus line -- DEFAULT IS NONE
726	write(,) "Titus crocus line ?"
727	tituscap=.false. ! default value
728	call getin("tituscap",tituscap)
729	write(,) "tituscap",tituscap
730
731	! Chemistry:
732	write(,) "photochemistry: include chemical species"
733	photochem=.false. ! default value
734	call getin("photochem",photochem)
735	write(,) " photochem = ",photochem
736
737	write(,) "Compute chemistry (if photochem is .true.)",
738	& "every ichemistry physics step (default: ichemistry=1)"
739	ichemistry=1
740	call getin("ichemistry",ichemistry)
741	write(,) " ichemistry = ",ichemistry
742
743
744	! SCATTERERS
745	write(,) "how many scatterers?"
746	naerkind=1 ! default value
747	call getin("naerkind",naerkind)
748	write(,)" naerkind = ",naerkind
749
750	! Test of incompatibility
751	c Logical tests for radiatively active water-ice clouds:
752	IF ( (activice.AND.(.NOT.water)).OR.
753	& (activice.AND.(naerkind.LT.2)) ) THEN
754	WRITE(,) 'If activice is TRUE, water has to be set'
755	WRITE(,) 'to TRUE, and "naerkind" must be at least'
756	WRITE(,) 'equal to 2.'
757	CALL ABORT
758	ENDIF
759
760	!------------------------------------------
761	!------------------------------------------
762	! once naerkind is known allocate arrays
763	! -- we do it here and not in phys_var_init
764	! -- because we need to know naerkind first
765	CALL ini_scatterers(ngrid,nlayer)
766	!------------------------------------------
767	!------------------------------------------
768
769
770	c Please name the different scatterers here ----------------
771	name_iaer(1) = "dust_conrath" !! default choice is good old Conrath profile
772	IF (doubleq.AND.active) name_iaer(1) = "dust_doubleq" !! two-moment scheme
773	if (nq.gt.1) then
774	! trick to avoid problems compiling with 1 tracer
775	! and picky compilers who know name_iaer(2) is out of bounds
776	j=2
777	IF (rdstorm.AND..NOT.activice.AND..NOT.slpwind) name_iaer(2) =
778	& "stormdust_doubleq" !! storm dust two-moment scheme
779	IF (rdstorm.AND.water.AND.activice.AND..NOT.slpwind)
780	& name_iaer(3) = "stormdust_doubleq"
781	IF (slpwind.AND..NOT.activice.AND..NOT.rdstorm) name_iaer(2) =
782	& "topdust_doubleq" !! storm dust two-moment scheme
783	IF (slpwind.AND.water.AND.activice.AND..NOT.rdstorm)
784	& name_iaer(3) = "topdust_doubleq"
785	IF (rdstorm.AND.slpwind.AND..NOT.activice) THEN
786	name_iaer(2) = "stormdust_doubleq"
787	name_iaer(3) = "topdust_doubleq"
788	ENDIF
789	IF (rdstorm.AND.slpwind.AND.water.AND.activice) THEN
790	name_iaer(3) = "stormdust_doubleq"
791	name_iaer(4) = "topdust_doubleq"
792	ENDIF
793	IF (water.AND.activice) name_iaer(j) = "h2o_ice" !! radiatively-active clouds
794	IF (submicron.AND.active) name_iaer(j) = "dust_submicron" !! JBM experimental stuff
795	endif ! of if (nq.gt.1)
796
797	c ----------------------------------------------------------
798
799	! THERMOSPHERE
800
801	write(,) "call thermosphere ?"
802	callthermos=.false. ! default value
803	call getin("callthermos",callthermos)
804	write(,) " callthermos = ",callthermos
805
806
807	write(,) " water included without cycle ",
808	& "(only if water=.false.)"
809	thermoswater=.false. ! default value
810	call getin("thermoswater",thermoswater)
811	write(,) " thermoswater = ",thermoswater
812
813	write(,) "call thermal conduction ?",
814	& " (only if callthermos=.true.)"
815	callconduct=.false. ! default value
816	call getin("callconduct",callconduct)
817	write(,) " callconduct = ",callconduct
818
819	write(,) "call EUV heating ?",
820	& " (only if callthermos=.true.)"
821	calleuv=.false. ! default value
822	call getin("calleuv",calleuv)
823	write(,) " calleuv = ",calleuv
824
825	write(,) "call molecular viscosity ?",
826	& " (only if callthermos=.true.)"
827	callmolvis=.false. ! default value
828	call getin("callmolvis",callmolvis)
829	write(,) " callmolvis = ",callmolvis
830
831	write(,) "call molecular diffusion ?",
832	& " (only if callthermos=.true.)"
833	callmoldiff=.false. ! default value
834	call getin("callmoldiff",callmoldiff)
835	write(,) " callmoldiff = ",callmoldiff
836
837
838	write(,) "call thermospheric photochemistry ?",
839	& " (only if callthermos=.true.)"
840	thermochem=.false. ! default value
841	call getin("thermochem",thermochem)
842	write(,) " thermochem = ",thermochem
843
844	write(,) "Method to include solar variability"
845	write(,) "0-> fixed value of E10.7 (fixed_euv_value); ",
846	& "1-> daily evolution of E10.7 (for given solvaryear)"
847	solvarmod=1
848	call getin("solvarmod",solvarmod)
849	write(,) " solvarmod = ",solvarmod
850
851	write(,) "Fixed euv (for solvarmod==0) 10.7 value?"
852	write(,) " (min=80 , ave=140, max=320)"
853	fixed_euv_value=140 ! default value
854	call getin("fixed_euv_value",fixed_euv_value)
855	write(,) " fixed_euv_value = ",fixed_euv_value
856
857	write(,) "Solar variability as observed for MY: "
858	write(,) "Only if solvarmod=1"
859	solvaryear=24
860	call getin("solvaryear",solvaryear)
861	write(,) " solvaryear = ",solvaryear
862
863	write(,) "UV heating efficiency:",
864	& "measured values between 0.19 and 0.23 (Fox et al. 1996)",
865	& "lower values may be used to compensate low 15 um cooling"
866	euveff=0.21 !default value
867	call getin("euveff",euveff)
868	write(,) " euveff = ", euveff
869
870
871	if (.not.callthermos) then
872	if (thermoswater) then
873	print*,'if thermoswater is set, callthermos must be true'
874	stop
875	endif
876	if (callconduct) then
877	print*,'if callconduct is set, callthermos must be true'
878	stop
879	endif
880	if (calleuv) then
881	print*,'if calleuv is set, callthermos must be true'
882	stop
883	endif
884	if (callmolvis) then
885	print*,'if callmolvis is set, callthermos must be true'
886	stop
887	endif
888	if (callmoldiff) then
889	print*,'if callmoldiff is set, callthermos must be true'
890	stop
891	endif
892	if (thermochem) then
893	print*,'if thermochem is set, callthermos must be true'
894	stop
895	endif
896	endif
897
898	! Test of incompatibility:
899	! if photochem is used, then water should be used too
900
901	if (photochem.and..not.water) then
902	print*,'if photochem is used, water should be used too'
903	stop
904	endif
905
906	! if callthermos is used, then thermoswater should be used too
907	! (if water not used already)
908
909	if (callthermos .and. .not.water) then
910	if (callthermos .and. .not.thermoswater) then
911	print*,'if callthermos is used, water or thermoswater
912	& should be used too'
913	stop
914	endif
915	endif
916
917	PRINT*,'--------------------------------------------'
918	PRINT*
919	PRINT*
920	ELSE
921	write(,)
922	write(,) 'Cannot read file callphys.def. Is it here ?'
923	stop
924	ENDIF
925
926	8000 FORMAT(t5,a12,l8)
927	8001 FORMAT(t5,a12,i8)
928
929	PRINT*
930	PRINT*,'conf_phys: daysec',daysec
931	PRINT*
932	PRINT*,'conf_phys: The radiative transfer is computed:'
933	PRINT*,' each ',iradia,' physical time-step'
934	PRINT,' or each ',iradiadtphys,' seconds'
935	PRINT*
936	! --------------------------------------------------------------
937	! Managing the Longwave radiative transfer
938	! --------------------------------------------------------------
939
940	! In most cases, the run just use the following values :
941	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
942	callemis=.true.
943	! ilwd=10*int(daysec/dtphys) ! bug before 22/10/01
944	ilwd=1
945	ilwn=1 !2
946	ilwb=1 !2
947	linear=.true.
948	ncouche=3
949	alphan=0.4
950	semi=0
951
952	! BUT people working hard on the LW may want to read them in 'radia.def'
953	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
954
955	OPEN(99,file='radia.def',status='old',form='formatted'
956	. ,iostat=ierr)
957	IF(ierr.EQ.0) THEN
958	write(,) 'conf_phys: Reading radia.def !!!'
959	READ(99,fmt='(a)') ch1
960	READ(99,*) callemis
961	WRITE(*,8000) ch1,callemis
962
963	READ(99,fmt='(a)') ch1
964	READ(99,*) iradia
965	WRITE(*,8001) ch1,iradia
966
967	READ(99,fmt='(a)') ch1
968	READ(99,*) ilwd
969	WRITE(*,8001) ch1,ilwd
970
971	READ(99,fmt='(a)') ch1
972	READ(99,*) ilwn
973	WRITE(*,8001) ch1,ilwn
974
975	READ(99,fmt='(a)') ch1
976	READ(99,*) linear
977	WRITE(*,8000) ch1,linear
978
979	READ(99,fmt='(a)') ch1
980	READ(99,*) ncouche
981	WRITE(*,8001) ch1,ncouche
982
983	READ(99,fmt='(a)') ch1
984	READ(99,*) alphan
985	WRITE(,) ch1,alphan
986
987	READ(99,fmt='(a)') ch1
988	READ(99,*) ilwb
989	WRITE(*,8001) ch1,ilwb
990
991
992	READ(99,fmt='(a)') ch1
993	READ(99,'(l1)') callg2d
994	WRITE(*,8000) ch1,callg2d
995
996	READ(99,fmt='(a)') ch1
997	READ(99,*) semi
998	WRITE(,) ch1,semi
999	end if
1000	CLOSE(99)
1001
1002	END

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

Download in other formats:

Original Format