Context Navigation

conf_phys.F @ 2179

Last change on this file since 2179 was 2179, checked in by mvals, 5 years ago
Mars GCM: Add the sublimation/condensation latent heat release from the surface water ice in vdifc_mod.F. It can be activated or deactivated with the new flag "latentheat" (for now latentheat=false by default). MV
File size: 34.0 KB

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

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

Download in other formats:

Original Format