Context Navigation

conf_phys.F @ 2124

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

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

Download in other formats:

Original Format