Context Navigation

conf_phys.F @ 2156

Last change on this file since 2156 was 2149, checked in by emillour, 5 years ago
Mars GCM: Add F.Lott's non-orographic GW parametrization. Disabled by default for now, activated by setting calllott_nonoro=.true. in callphys.def GG+EM
File size: 32.9 KB

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

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

Download in other formats:

Original Format