Context Navigation

conf_phys.F @ 2282

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

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

Download in other formats:

Original Format