Context Navigation

conf_phys.F @ 4054

Last change on this file since 4054 was 4054, checked in by tbertrand, 11 days ago
Mars PCM: Adding coagulation of dust particles (cf Bertrand et al., 2022, Icarus) TB
File size: 52.9 KB

Line
1	MODULE conf_phys_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	SUBROUTINE conf_phys(ngrid,nlayer,nq)
8
9	!=======================================================================
10	!
11	! purpose:
12	! -------
13	!
14	! Initialisation for the physical parametrisations
15	! flags (i.e. run-time options) of the Mars PCM
16	!-----------------------------------------------------------------------
17
18	USE ioipsl_getin_p_mod, ONLY : getin_p
19	use tracer_mod, only : nuice_sed, ccn_factor, nuiceco2_sed,
20	& nuice_ref,nuiceco2_ref
21	use surfdat_h, only: albedo_h2o_cap,albedo_h2o_frost,
22	& albedo_perennialco2,
23	& frost_albedo_threshold, inert_h2o_ice,
24	& frost_metam_threshold,old_wsublimation_scheme
25	use time_phylmdz_mod, only: steps_per_sol,outputs_per_sol,iphysiq,
26	& ecritstart,daysec,dtphys,slow_diagfi
27	use dimradmars_mod, only: naerkind, name_iaer,
28	& ini_scatterers,tauvis
29	use datafile_mod, only: datadir
30	use wstats_mod, only: callstats
31	use writediagsoil_mod, only: diagsoil
32	use calchim_mod, only: ichemistry
33	use co2condens_mod, only: scavco2cond
34	use dust_param_mod, only: dustbin, doubleq, submicron, active,
35	& lifting, freedust, callddevil,
36	& dustscaling_mode,
37	& reff_driven_IRtoVIS_scenario
38	use aeropacity_mod, only: iddist, topdustref
39	#ifndef MESOSCALE
40	USE mod_phys_lmdz_transfert_para, ONLY: bcast
41	#endif
42	USE paleoclimate_mod,ONLY: paleoclimate,
43	& lag_layer, include_waterbuoyancy
44	use microphys_h, only: mteta
45	use comsoil_h, only: adsorption_soil, choice_ads,ads_const_D,
46	& ads_massive_ice
47	use nonoro_gwd_mix_mod, only: calljliu_gwimix
48	use callkeys_mod, only: startphy_file, activice, activeco2ice,
49	& alphan, calladj, callatke, callcond,
50	& callconduct, calldifv, callemis,
51	& calleofdump, calleuv, callg2d, calllott,
52	& calllott_nonoro, callmoldiff, callmolvis,
53	& callnirco2, callnlte, callrad, callrichsl,
54	& callslope, callsoil, calltherm, caps,
55	& callthermos, callyamada4, CLFvarying,
56	& CLFvaryingCO2, CLFfixval,
57	& cloud_adapt_ts, co2clouds,
58	& co2useh2o, coeff_detrainment, diurnal,
59	& coeff_injection, dustinjection,
60	& cst_cap_albedo, dustiropacity, dustrefir,
61	& euveff, fixed_euv_value, hdo, hdofrac,
62	& iaervar, ilwb, ilwd, ilwn, iradia,
63	& linear, latentheat_surfwater, microphys,
64	& meteo_flux, ncouche, nircorr, nltemodel,
65	& photochem, poreice_tifeedback, rayleigh,
66	& rdstorm, refill_watercap, satindexco2,
67	& scavenging, season, sedimentation, semi,
68	& solvarmod, solvaryear, spantCO2, supersat,
69	& surfaceice_tifeedback, temp_dependent_m,
70	& tf_injection, ti_injection, thermochem,
71	& thermoswater, tituscap, tke_heat_flux,
72	& topflows, water, coagulation, kernel_b,
73	& kernel_g,kernel_de,kernel_ti,coal_kg,
74	& fullcoag
75
76	IMPLICIT NONE
77
78	INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
79	INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
80	INTEGER,INTENT(IN) :: nq ! number of tracers
81
82	INTEGER ierr,j
83	character(len=20),parameter :: modname="conf_phys"
84
85	real :: ecritphy ! to check that this obsolete flag is no longer used...
86	real :: albedo_perennialco2_north, albedo_perennialco2_south
87
88	CHARACTER ch1*12
89	#ifndef MESOSCALE
90	! read in some parameters from "run.def" for physics,
91	! or shared between dynamics and physics.
92	ecritphy=-666 ! dummy default value
93	call getin_p("ecritphy",ecritphy) ! frequency of outputs in physics,
94	! in dynamical steps
95	if (ecritphy/=-666) then
96	call abort_physic(modname,
97	& "Error: parameter ecritphy is obsolete! Remove it! "//
98	& "And use outputs_per_sol instead",1)
99	endif
100
101	outputs_per_sol=4 ! default value, 4 outputs per sol
102	call getin_p("outputs_per_sol",outputs_per_sol)
103	! check that indeed it is possible for given physics time step
104	if (modulo(steps_per_sol,outputs_per_sol)/=0) then
105	write(,) "Error: outputs_per_sol = ",outputs_per_sol
106	write(,) " is not a divisor of number of steps per sol = ",
107	& steps_per_sol
108	call abort_physic(modname,
109	& "Error: inadequate value for outputs_per_sol",1)
110	endif
111
112	slow_diagfi =.false. ! default value, open/close diagfi once (1D only)
113	call getin_p("slow_diagfi",slow_diagfi) ! if true, open/close the diagfi.nc
114	! at output frequency instead of once
115	! (slower but required for debug)
116	write(,) "Develop: slow_diagfi is ", slow_diagfi
117
118	iphysiq=20 ! default value
119	call getin_p("iphysiq",iphysiq) ! call physics every iphysiq dyn step
120	ecritstart=0 ! default value
121	call getin_p("ecritstart",ecritstart) ! write a restart every ecristart steps
122	#endif
123
124	! --------------------------------------------------------------
125	! Reading the "callphys.def" file controlling some key options
126	! --------------------------------------------------------------
127
128	!$OMP MASTER
129	! check that 'callphys.def' file is around
130	OPEN(99,file='callphys.def',status='old',form='formatted'
131	& ,iostat=ierr)
132	CLOSE(99)
133	!$OMP END MASTER
134	#ifndef MESOSCALE
135	call bcast(ierr)
136	#endif
137	! ierr=0
138
139	IF(ierr.EQ.0) THEN
140	PRINT*
141	PRINT*
142	PRINT*,'--------------------------------------------'
143	PRINT*,' conf_phys: Parameters for the physics (callphys.def)'
144	PRINT*,'--------------------------------------------'
145
146	write(,) "Directory where external input files are:"
147	! default path is set in datafile_mod
148	call getin_p("datadir",datadir)
149	write(,) " datadir = ",trim(datadir)
150
151	write(,) "Initialize physics with startfi.nc file ?"
152	startphy_file=.true.
153	call getin_p("startphy_file",startphy_file)
154	write(,) "startphy_file", startphy_file
155
156	write(,) "Diurnal cycle ?"
157	write(,) "(if diurnal=False, diurnal averaged solar heating)"
158	diurnal=.true. ! default value
159	call getin_p("diurnal",diurnal)
160	write(,) " diurnal = ",diurnal
161
162	write(,) "Seasonal cycle ?"
163	write(,) "(if season=False, Ls stays constant, to value ",
164	& "set in 'start'"
165	season=.true. ! default value
166	call getin_p("season",season)
167	write(,) " season = ",season
168
169	write(,) "Save statistics in file stats.nc ?"
170	#ifdef MESOSCALE
171	callstats=.false. ! default value
172	#else
173	callstats=.true. ! default value
174	#endif
175	call getin_p("callstats",callstats)
176	write(,) " callstats = ",callstats
177
178	write(,) "Write sub-surface fields in file diagsoil.nc ?"
179	diagsoil=.false. ! default value
180	call getin_p("diagsoil",diagsoil)
181	write(,) " diagsoil = ",diagsoil
182
183	write(,) "Save EOF profiles in file 'profiles' for ",
184	& "Climate Database?"
185	calleofdump=.false. ! default value
186	call getin_p("calleofdump",calleofdump)
187	write(,) " calleofdump = ",calleofdump
188
189	write(,) "Dust scenario: 1=constant dust (read from startfi",
190	& " or set as tauvis); 2=Viking scenario; =3 MGS scenario,",
191	& "=6 cold (low dust) scenario; =7 warm (high dust) scenario ",
192	& "=24,25 ... 30 :Mars Year 24, ... or 30 from TES assimilation"
193	iaervar=3 ! default value
194	call getin_p("iaervar",iaervar)
195	write(,) " iaervar = ",iaervar
196
197	write(,) "Reference (visible) dust opacity at 610 Pa ",
198	& "(matters only if iaervar=1)"
199	! NB: default value of tauvis is set/read in startfi.nc file
200	call getin_p("tauvis",tauvis)
201	write(,) " tauvis = ",tauvis
202
203	write(,) "Dust vertical distribution:"
204	write(,) "(=1 top set by topdustref parameter;",
205	& " =2 Viking scenario; =3 MGS scenario)"
206	iddist=3 ! default value
207	call getin_p("iddist",iddist)
208	write(,) " iddist = ",iddist
209
210	write(,) "Dust top altitude (km). (Matters only if iddist=1)"
211	topdustref= 90.0 ! default value
212	call getin_p("topdustref",topdustref)
213	write(,) " topdustref = ",topdustref
214
215	write(,) "Prescribed surface thermal flux (H/(rho*cp),K m/s)"
216	tke_heat_flux=0. ! default value
217	call getin_p("tke_heat_flux",tke_heat_flux)
218	write(,) " tke_heat_flux = ",tke_heat_flux
219	write(,) " 0 means the usual schemes are computing"
220
221	write(,) "call radiative transfer ?"
222	callrad=.true. ! default value
223	call getin_p("callrad",callrad)
224	write(,) " callrad = ",callrad
225
226	write(,) "call slope insolation scheme ?",
227	& "(matters only if callrad=T)"
228	#ifdef MESOSCALE
229	callslope=.true. ! default value
230	#else
231	callslope=.false. ! default value (not supported yet)
232	#endif
233	call getin_p("callslope",callslope)
234	write(,) " callslope = ",callslope
235
236	write(,) "call NLTE radiative schemes ?",
237	& "(matters only if callrad=T)"
238	callnlte=.false. ! default value
239	call getin_p("callnlte",callnlte)
240	write(,) " callnlte = ",callnlte
241
242	nltemodel=0 !default value
243	write(,) "NLTE model?"
244	write(,) "0 -> old model, static O"
245	write(,) "1 -> old model, dynamic O"
246	write(,) "2 -> new model"
247	write(,) "(matters only if callnlte=T)"
248	call getin_p("nltemodel",nltemodel)
249	write(,) " nltemodel = ",nltemodel
250
251	write(,) "call CO2 NIR absorption ?",
252	& "(matters only if callrad=T)"
253	callnirco2=.false. ! default value
254	call getin_p("callnirco2",callnirco2)
255	write(,) " callnirco2 = ",callnirco2
256
257	write(,) "New NIR NLTE correction ?",
258	$ "0-> old model (no correction)",
259	$ "1-> new correction",
260	$ "(matters only if callnirco2=T)"
261
262	nircorr=0 !default value. this is OK below 60 km.
263	! nircorr=1 should be prefered, but requires that co2
264	! and o tracers be available
265	call getin_p("nircorr",nircorr)
266	write(,) " nircorr = ",nircorr
267
268	write(,) "call turbulent vertical diffusion ?"
269	calldifv=.true. ! default value
270	call getin_p("calldifv",calldifv)
271	write(,) " calldifv = ",calldifv
272
273	write(,) "call thermals ?"
274	calltherm=.false. ! default value
275	call getin_p("calltherm",calltherm)
276	write(,) " calltherm = ",calltherm
277
278	write(,) "call convective adjustment ?"
279	calladj=.true. ! default value
280	call getin_p("calladj",calladj)
281	write(,) " calladj = ",calladj
282
283	if (calltherm .and. calladj) then
284	print*,'!!! PLEASE NOTE !!!'
285	print*,'convective adjustment is on'
286	print*,'but since thermal plume model is on'
287	print*,'convadj is only activated above the PBL'
288	endif
289
290	write(,) "used latest version of yamada scheme?"
291	callyamada4=.true. ! default value
292	call getin_p("callyamada4",callyamada4)
293	write(,) " callyamada4 = ",callyamada4
294
295	write(,) "used latest version of ATKE scheme?"
296	callatke =.false. ! default value
297	call getin_p("callatke",callatke)
298	write(,) " callatke = ",callatke
299
300	if(callyamada4.and.callatke) then
301	print*,"You can't use both yamada and atke"
302	print*,"Choose either Yamada or ATKE"
303	call abort_physic(modname,
304	& "Can't use both yamada and ATKE",1)
305	endif
306
307	if (calltherm .and. .not.callyamada4) then
308	print*,'!!!! WARNING WARNING WARNING !!!!'
309	print*,'if calltherm=T we strongly advise that '
310	print*,'you set the flag callyamada4 to T '
311	print*,'!!!! WARNING WARNING WARNING !!!!'
312	endif
313
314	write(,) "call Richardson-based surface layer ?"
315	callrichsl=.false. ! default value
316	call getin_p("callrichsl",callrichsl)
317	write(,) " callrichsl = ",callrichsl
318
319	if (calltherm .and. .not.callrichsl) then
320	print*,'WARNING WARNING WARNING'
321	print*,'if calltherm=T we strongly advise that '
322	print*,'you use the new surface layer scheme '
323	print*,'by setting callrichsl=T '
324	endif
325
326	if (calladj .and. callrichsl .and. (.not. calltherm)) then
327	print*,'You should not be calling the convective adjustment
328	& scheme with the Richardson surface-layer and without the thermals
329	&. This approach is not
330	& physically consistent and can lead to unrealistic friction
331	& values.'
332	print*,'If you want to use the Ri. surface-layer, either
333	& activate thermals OR de-activate the convective adjustment.'
334	call abort_physic(modname,
335	& "Richardson layer must be used with thermals",1)
336	endif
337
338	write(,) "call CO2 condensation ?"
339	callcond=.true. ! default value
340	call getin_p("callcond",callcond)
341	write(,) " callcond = ",callcond
342
343	write(,)"call thermal conduction in the soil ?"
344	callsoil=.true. ! default value
345	call getin_p("callsoil",callsoil)
346	write(,) " callsoil = ",callsoil
347
348
349	write(,)"call Lott's gravity wave/subgrid topography ",
350	& "scheme ?"
351	calllott=.true. ! default value
352	call getin_p("calllott",calllott)
353	write(,)" calllott = ",calllott
354
355	write(,)"call Lott's non-oro GWs parameterisation ",
356	& "scheme ?"
357	calllott_nonoro=.false. ! default value
358	call getin_p("calllott_nonoro",calllott_nonoro)
359	write(,)" calllott_nonoro = ",calllott_nonoro
360
361	write(,)"call jliu's nonoro GWs-induced mixing ",
362	& "scheme ?"
363	calljliu_gwimix=.false. ! default value
364	call getin_p("calljliu_gwimix",calljliu_gwimix)
365	write(,)" calljliu_gwimix = ",calljliu_gwimix
366
367
368	! rocket dust storm injection scheme
369	write(,)"call rocket dust storm parametrization"
370	rdstorm=.false. ! default value
371	call getin_p("rdstorm",rdstorm)
372	write(,)" rdstorm = ",rdstorm
373	! rocket dust storm detrainment coefficient
374	coeff_detrainment=0.02 ! default value
375	call getin_p("coeff_detrainment",coeff_detrainment)
376	write(,)" coeff_detrainment = ",coeff_detrainment
377
378	! entrainment by mountain top dust flows scheme
379	write(,)"call mountain top dust flows parametrization"
380	topflows=.false. ! default value
381	call getin_p("topflows",topflows)
382	write(,)" topflows = ",topflows
383
384	! dust particle coagulation
385	write(,)"call coagulation of dust"
386	coagulation=.false. ! default value
387	call getin_p("coagulation",coagulation)
388	write(,)" coagulation = ",coagulation
389	fullcoag=.false. ! default value
390	kernel_b=.true. ! default value (if coagulation=true)
391	kernel_g=.false. ! default value
392	kernel_de=.false. ! default value
393	kernel_ti=.false. ! default value
394	coal_kg=0 ! default value
395	if (coagulation) then
396	write(,)" coagulation fullcoag and mode= ",fullcoag,coal_kg
397	write(,)" coagulation kernels= ",kernel_b,kernel_g,
398	& kernel_de,kernel_ti
399	endif
400
401	! latent heat release from ground water ice sublimation/condensation
402	write(,)"latent heat release during sublimation",
403	& " /condensation of ground water ice"
404	latentheat_surfwater=.true. ! default value
405	call getin_p("latentheat_surfwater",latentheat_surfwater)
406	write(,)" latentheat_surfwater = ",latentheat_surfwater
407
408	write(,)"rad.transfer is computed every iradia",
409	& " physical timestep"
410	iradia=1 ! default value
411	call getin_p("iradia",iradia)
412	write(,)" iradia = ",iradia
413
414
415	write(,)"Output of the exchange coefficient mattrix ?",
416	& "(for diagnostics only)"
417	callg2d=.false. ! default value
418	call getin_p("callg2d",callg2d)
419	write(,)" callg2d = ",callg2d
420
421	write(,)"Rayleigh scattering : (should be .false. for now)"
422	rayleigh=.false.
423	call getin_p("rayleigh",rayleigh)
424	write(,)" rayleigh = ",rayleigh
425
426	! PALEOCLIMATE
427
428	write(,)"Using lag layer??"
429	lag_layer=.false.
430	call getin_p("lag_layer",lag_layer)
431	write(,) "lag_layer = ", lag_layer
432
433	write(,)"Is it paleoclimate run?"
434	paleoclimate=.false. ! default value
435	call getin_p("paleoclimate",paleoclimate)
436	write(,)"paleoclimate = ",paleoclimate
437
438	write(,)"Albedo for perennial CO2 ice (northern hemisphere)?"
439	albedo_perennialco2_north = 0.6 ! default value
440	call getin_p("albedo_perennialco2_north",albedo_perennialco2_
441	& north)
442	albedo_perennialco2(1) = albedo_perennialco2_north
443	write(,)"albedo_perennialco2(1) = ",albedo_perennialco2(1)
444
445	write(,)"Albedo for perennial CO2 ice (southern hemisphere)?"
446	albedo_perennialco2_south = 0.85 ! default value
447	call getin_p("albedo_perennialco2_south",albedo_perennialco2_
448	& south)
449	albedo_perennialco2(2) = albedo_perennialco2_south
450	write(,)"albedo_perennialco2(2) = ",albedo_perennialco2(2)
451
452	write(,)"Include water buoyancy effect??"
453	include_waterbuoyancy=.false.
454	call getin_p("include_waterbuoyancy",include_waterbuoyancy)
455	write(,) "include_waterbuoyancy = ", include_waterbuoyancy
456
457	if (include_waterbuoyancy.and.(.not.paleoclimate)) then
458	print*,'You should not include the water buoyancy effect
459	& on current climate (model not tunned with this effect'
460	call abort_physic(modname,
461	& "include_waterbuoyancy must be used with paleoclimate",1)
462	endif
463
464	! TRACERS:
465
466	! dustbin
467	write(,)"Transported dust ? (if >0, use 'dustbin' dust bins)"
468	dustbin=0 ! default value
469	call getin_p("dustbin",dustbin)
470	write(,)" dustbin = ",dustbin
471	! active
472	write(,)"Radiatively active dust ? (matters if dustbin>0)"
473	active=.false. ! default value
474	call getin_p("active",active)
475	write(,)" active = ",active
476
477	! Test of incompatibility:
478	! if active is used, then dustbin should be > 0
479
480	if (active.and.(dustbin.lt.1)) then
481	print*,'if active is used, then dustbin should > 0'
482	call abort_physic(modname,
483	& "active option requires dustbin > 0",1)
484	endif
485	! doubleq
486	write(,)"use mass and number mixing ratios to predict",
487	& " dust size ?"
488	doubleq=.false. ! default value
489	call getin_p("doubleq",doubleq)
490	write(,)" doubleq = ",doubleq
491	! submicron
492	submicron=.false. ! default value
493	call getin_p("submicron",submicron)
494	write(,)" submicron = ",submicron
495
496	! Test of incompatibility:
497	! if doubleq is used, then dustbin should be 2
498
499	if (doubleq.and.(dustbin.ne.2)) then
500	print*,'if doubleq is used, then dustbin should be 2'
501	call abort_physic(modname,
502	& "doubleq option requires dustbin = 2",1)
503	endif
504	if (doubleq.and.submicron.and.(nq.LT.3)) then
505	print*,'If doubleq is used with a submicron tracer,'
506	print*,' then the number of tracers has to be'
507	print*,' larger than 3.'
508	call abort_physic(modname,
509	& "submicron option requires dustbin > 2",1)
510	endif
511
512	! lifting
513	write(,)"dust lifted by GCM surface winds ?"
514	lifting=.false. ! default value
515	call getin_p("lifting",lifting)
516	write(,)" lifting = ",lifting
517
518	! Test of incompatibility:
519	! if lifting is used, then dustbin should be > 0
520
521	if (lifting.and.(dustbin.lt.1)) then
522	print*,'if lifting is used, then dustbin should > 0'
523	call abort_physic(modname,
524	& "lifting option requires dustbin > 0",1)
525	endif
526
527	! dust injection scheme
528	dustinjection=0 ! default: no injection scheme
529	call getin_p("dustinjection",dustinjection)
530	write(,)" dustinjection = ",dustinjection
531	! dust injection scheme coefficient
532	coeff_injection=0.25 ! default value
533	call getin_p("coeff_injection",coeff_injection)
534	write(,)" coeff_in,jection = ",coeff_injection
535	! timing for dust injection
536	ti_injection=0. ! default value
537	tf_injection=24. ! default value
538	call getin_p("ti_injection",ti_injection)
539	write(,)" ti_injection = ",ti_injection
540	call getin_p("tf_injection",tf_injection)
541	write(,)" tf_injection = ",tf_injection
542
543	! free evolving dust
544	! freedust=true just says that there is no lifting and no dust opacity scaling.
545	write(,)"dust lifted by GCM surface winds ?"
546	freedust=.false. ! default value
547	call getin_p("freedust",freedust)
548	write(,)" freedust = ",freedust
549	if (freedust.and..not.doubleq) then
550	print*,'freedust should be used with doubleq !'
551	call abort_physic(modname,
552	& "freedust option requires doubleq",1)
553	endif
554
555	! dust rescaling mode (if any)
556	! =0 --> freedust, tauscaling=1
557	! =1 --> GCM5.3-like, tauscaling computed wrt tau_pref_scenario
558	! =2 --> tauscaling=1 but dust_rad_adjust (rescaling for the radiative transfer only)
559	if (freedust) then
560	dustscaling_mode=0
561	else
562	dustscaling_mode=1 ! GCMv5.3 style
563	endif
564	call getin_p("dustscaling_mode",dustscaling_mode)
565	write(,) "dustscaling_mode=",dustscaling_mode
566	if ((dustscaling_mode.eq.1).and.freedust) then
567	print*,'freedust and dustscaling_mode=1 are incompatible !'
568	call abort_physic(modname,
569	& "freedust and dustscaling_mode=1 are incompatible",1)
570	endif
571
572	#ifndef MESOSCALE
573	! this test is valid in GCM case
574	! ... not in mesoscale case, for we want to activate mesoscale lifting
575	if (freedust.and.dustinjection.eq.0)then
576	if(lifting) then
577	print*,'if freedust is used and dustinjection = 0,
578	& then lifting should not be used'
579	call abort_physic(modname,
580	& "freedust option with dustinjection = 0"//
581	& " requires lifting to be false",1)
582	endif
583	endif
584	#endif
585	if (dustinjection.eq.1)then
586	if(.not.lifting) then
587	print*,"if dustinjection=1, then lifting should be true"
588	call abort_physic(modname,
589	& "dustinjection=1 requires lifting",1)
590	endif
591	if(.not.freedust) then
592	print*,"if dustinjection=1, then freedust should be true"
593	call abort_physic(modname,
594	& "dustinjection=1 requires freedust",1)
595	endif
596	endif
597	! rocket dust storm and entrainment by top flows
598	! Test of incompatibility:
599	! if rdstorm or topflows is used, then doubleq should be true
600	if ((rdstorm.or.topflows).and..not.doubleq) then
601	print*,'if rdstorm or topflows is used, then doubleq
602	& should be used !'
603	call abort_physic(modname,
604	& "rdstorm or topflows requires doubleq",1)
605	endif
606	if ((rdstorm.or.topflows).and..not.active) then
607	print*,'if rdstorm or topflows is used, then active
608	& should be used !'
609	call abort_physic(modname,
610	& "rdstorm or topflows requires active",1)
611	endif
612	if (rdstorm.and..not.lifting) then
613	print*,'if rdstorm is used, then lifting
614	& should be used !'
615	call abort_physic(modname,
616	& "rdstorm requires lifting",1)
617	endif
618	if ((rdstorm.or.topflows).and..not.freedust) then
619	print*,'if rdstorm or topflows is used, then freedust
620	& should be used !'
621	call abort_physic(modname,
622	& "rdstorm or topflows requires freedust",1)
623	endif
624	if (rdstorm.and.(dustinjection.eq.0)) then
625	print*,'if rdstorm is used, then dustinjection
626	& should be used !'
627	call abort_physic(modname,
628	& "rdstorm requires dustinjection",1)
629	endif
630	! Dust IR opacity
631	write(,)" Wavelength for infrared opacity of dust ?"
632	write(,)" Choices are:"
633	write(,)" tes --- > 9.3 microns [default]"
634	write(,)" mcs --- > 21.6 microns"
635	!
636	! WARNING WARNING WARNING WARNING WARNING WARNING
637	!
638	! BEFORE ADDING A NEW VALUE, BE SURE THAT THE
639	! CORRESPONDING WAVELENGTH IS IN THE LOOKUP TABLE,
640	! OR AT LEAST NO TO FAR, TO AVOID FALLACIOUS INTERPOLATIONS.
641	!
642	dustiropacity="tes" !default value
643	call getin_p("dustiropacity",dustiropacity)
644	write(,)" dustiropacity = ",trim(dustiropacity)
645	select case (trim(dustiropacity))
646	case ("tes")
647	dustrefir = 9.3E-6
648	case ("mcs")
649	dustrefir = 21.6E-6
650	case default
651	write(,) trim(dustiropacity),
652	& " is not a valid option for dustiropacity"
653	call abort_physic(modname,
654	& "invalid dustiropacity option value",1)
655	end select
656	! Dust scenario IR to VIS conversion
657	write(,)"Use an IR to VIS conversion coefficient"
658	write(,)"for the dust scenario, that is dependent"
659	write(,)"on the GCM dust effective radius,"
660	write(,)"instead of a fixed 2.6 coefficient ?"
661	reff_driven_IRtoVIS_scenario=.false. !default value
662	call getin_p("reff_driven_IRtoVIS_scenario",
663	& reff_driven_IRtoVIS_scenario)
664	write(,)" reff_driven_IRtoVIS_scenario = ",
665	& reff_driven_IRtoVIS_scenario
666	! Test of incompatibility:
667	! if reff_driven_IRtoVIS_scenario=.true.,
668	! dustrefir must be 9.3E-6 = scenarios' wavelength
669	if (reff_driven_IRtoVIS_scenario .and.
670	& (dustrefir.ne.9.3E-6)) then
671	print*,'if reff_driven_IRtoVIS_scenario is used, then '//
672	& 'the GCM IR reference wavelength should be the one '//
673	& 'of the scenarios (dustiropacity=tes)'
674	call abort_physic(modname,
675	& "reff_driven_IRtoVIS_scenario requires tes wavelength",1)
676	endif
677
678	! callddevil
679	write(,)" dust lifted by dust devils ?"
680	callddevil=.false. !default value
681	call getin_p("callddevil",callddevil)
682	write(,)" callddevil = ",callddevil
683
684	! Test of incompatibility:
685	! if dustdevil is used, then dustbin should be > 0
686	if (callddevil.and.(dustbin.lt.1)) then
687	print*,'if dustdevil is used, then dustbin should > 0'
688	call abort_physic(modname,
689	& "callddevil requires dustbin > 0",1)
690	endif
691
692	! sedimentation
693	write(,) "Gravitationnal sedimentation ?"
694	sedimentation=.true. ! default value
695	call getin_p("sedimentation",sedimentation)
696	write(,) " sedimentation = ",sedimentation
697	! activice
698	write(,) "Radiatively active transported atmospheric ",
699	& "water ice ?"
700	activice=.false. ! default value
701	call getin_p("activice",activice)
702	write(,) " activice = ",activice
703	! water
704	write(,) "Compute water cycle ?"
705	water=.false. ! default value
706	call getin_p("water",water)
707	write(,) " water = ",water
708	! hdo
709	write(,) "Compute hdo cycle ?"
710	hdo=.false. ! default value
711	call getin_p("hdo",hdo)
712	write(,) " hdo = ",hdo
713
714	write(,) "Use fractionation for hdo?"
715	hdofrac=.true. ! default value
716	call getin_p("hdofrac",hdofrac)
717	write(,) " hdofrac = ",hdofrac
718
719	! Activeco2ice
720	write(,) "Radiatively active transported atmospheric ",
721	& "Co2 ice ?"
722	activeco2ice=.false. ! default value
723	call getin_p("activeco2ice",activeco2ice)
724	write(,) " activeco2ice = ",activeco2ice
725
726	! sub-grid cloud fraction: fixed
727	write(,) "Fixed cloud fraction?"
728	CLFfixval=1.0 ! default value
729	call getin_p("CLFfixval",CLFfixval)
730	write(,) "CLFfixval=",CLFfixval
731	! sub-grid cloud fraction: varying
732	write(,) "Use partial nebulosity?"
733	CLFvarying=.false. ! default value
734	call getin_p("CLFvarying",CLFvarying)
735	write(,)"CLFvarying=",CLFvarying
736
737	!CO2 clouds scheme?
738	write(,) "Compute CO2 clouds (implies microphysical scheme)?"
739	co2clouds=.false. ! default value
740	call getin_p("co2clouds",co2clouds)
741	write(,) " co2clouds = ",co2clouds
742	!Can water ice particles serve as CCN for CO2clouds
743	write(,) "Use water ice as CO2 clouds CCN ?"
744	co2useh2o=.false. ! default value
745	call getin_p("co2useh2o",co2useh2o)
746	write(,) " co2useh2o = ",co2useh2o
747	!Do we allow a supply of meteoritic paricles to serve as CO2 ice CCN?
748	write(,) "Supply meteoritic particle for CO2 clouds ?"
749	meteo_flux=.false. !Default value
750	call getin_p("meteo_flux",meteo_flux)
751	write(,) " meteo_flux = ",meteo_flux
752	!Do we allow a sub-grid temperature distribution for the CO2 microphysics
753	write(,) "sub-grid temperature distribution for CO2 clouds?"
754	CLFvaryingCO2=.false. !Default value
755	call getin_p("CLFvaryingCO2",CLFvaryingCO2)
756	write(,) " CLFvaryingCO2 = ",CLFvaryingCO2
757	!Amplitude of the sub-grid temperature distribution for the CO2 microphysics
758	write(,) "sub-grid temperature amplitude for CO2 clouds?"
759	spantCO2=0 !Default value
760	call getin_p("spantCO2",spantCO2)
761	write(,) " spantCO2 = ",spantCO2
762	!Do you want to filter the sub-grid T distribution by a Saturation index?
763	write(,) "filter sub-grid temperature by Saturation index?"
764	satindexco2=.true.
765	call getin_p("satindexco2",satindexco2)
766	write(,) " satindexco2 = ",satindexco2
767
768
769	! thermal inertia feedback
770	write(,) "Activate the thermal inertia feedback ?"
771	surfaceice_tifeedback=.false. ! default value
772	call getin_p("surfaceice_tifeedback",surfaceice_tifeedback)
773	write(,) " surfaceice_tifeedback = ",surfaceice_tifeedback
774
775	! Test of incompatibility:
776
777	if (surfaceice_tifeedback.and..not.water) then
778	print*,'if surfaceice_tifeedback is used,'
779	print*,'water should be used too'
780	call abort_physic(modname,
781	& "surfaceice_tifeedback requires water",1)
782	endif
783
784	if (surfaceice_tifeedback.and..not.callsoil) then
785	print*,'if surfaceice_tifeedback is used,'
786	print*,'callsoil should be used too'
787	call abort_physic(modname,
788	& "surfaceice_tifeedback requires callsoil",1)
789	endif
790
791	if (activice.and..not.water) then
792	print*,'if activice is used, water should be used too'
793	call abort_physic(modname,
794	& "activice requires water",1)
795	endif
796
797	if (hdo.and..not.water) then
798	print*,'if hdo is used, water should be used too'
799	call abort_physic(modname,
800	& "hd2 requires tracer",1)
801	endif
802
803
804	if (activeco2ice.and..not.co2clouds) then
805	print*,'if activeco2ice is used, co2clouds should be used too'
806	call abort_physic(modname,
807	& "activeco2ice requires co2clouds",1)
808	endif
809
810	! water ice clouds effective variance distribution for sedimentaion
811	write(,) "Sed effective variance for water ice clouds ?"
812	nuice_sed=0.45
813	call getin_p("nuice_sed",nuice_sed)
814	write(,) "water_param nueff Sedimentation:", nuice_sed
815
816	write(,) "Sed effective variance for CO2 clouds ?"
817	nuiceco2_sed=0.45
818	call getin_p("nuiceco2_sed",nuiceco2_sed)
819	write(,) "CO2 nueff Sedimentation:", nuiceco2_sed
820
821	write(,) "REF effective variance for CO2 clouds ?"
822	nuiceco2_ref=0.45
823	call getin_p("nuiceco2_ref",nuiceco2_ref)
824	write(,) "CO2 nueff Sedimentation:", nuiceco2_ref
825
826	write(,) "REF effective variance for water clouds ?"
827	nuice_ref=0.1
828	call getin_p("nuice_ref",nuice_ref)
829	write(,) "H2O nueff Sedimentation:", nuice_ref
830
831
832	! ccn factor if no scavenging
833	write(,) "water param CCN reduc. factor ?"
834	ccn_factor = 4.5
835	call getin_p("ccn_factor",ccn_factor)
836	write(,)" ccn_factor = ",ccn_factor
837	write(,)"Careful: only used when microphys=F, otherwise"
838	write(,)"the contact parameter is used instead;"
839
840	! microphys
841	write(,)"Microphysical scheme for water-ice clouds?"
842	microphys=.false. ! default value
843	call getin_p("microphys",microphys)
844	write(,)" microphys = ",microphys
845
846	! supersat
847	write(,)"Allow super-saturation of water vapor?"
848	supersat=.true. ! default value
849	call getin_p("supersat",supersat)
850	write(,)"supersat = ",supersat
851
852	! microphysical parameter contact
853	write(,) "water contact parameter ?"
854	mteta = 0.95 ! default value
855	temp_dependent_m = .false. ! default value
856	call getin_p("temp_dependent_m",temp_dependent_m)
857	if (temp_dependent_m) then !(JN 2023)
858	print*,'You have chosen a temperature-dependent water'
859	print*,'contact parameter ! From Maattanen et al. 2014'
860	else if (.not.temp_dependent_m) then
861	print*,'Water contact parameter is constant'
862	call getin_p("mteta",mteta)
863	write(,) "mteta = ", mteta
864	endif
865
866	! Adaptative timestep for cloud microphysics (JN 2023)
867	write(,)"Adaptative timestep for cloud",
868	& " microphysics ? (default is false)"
869	cloud_adapt_ts=.false. ! default value
870	call getin_p("cloud_adapt_ts",cloud_adapt_ts)
871	write(,)"cloud_adapt_ts= ",cloud_adapt_ts
872
873	! Test of incompatibility:
874	! If you want the adaptative timestep, you should use the
875	! temperature dependent contact parameter (otherwise the
876	! global water cycle is nuts !)
877	! However one can use the temperature dependent contact parameter
878	! without the adaptative timestep (MCD6.1 configuration)
879
880	if (cloud_adapt_ts.and..not.temp_dependent_m) then
881	print*,'Water cycle v6 : if cloud_adapt_ts is used'
882	print*,'then temp_dependent_m must be used!'
883	print*,'Otherwise the water cycle is unrealistic.'
884	call abort_physic(modname,
885	& "cloud_adapt_ts requires temp_dependent_m",1)
886	endif
887	! scavenging
888	write(,)"Dust scavenging by H2O/CO2 snowfall ?"
889	scavenging=.false. ! default value
890	call getin_p("scavenging",scavenging)
891	write(,)" scavenging = ",scavenging
892
893
894	! Test of incompatibility:
895	! if scavenging is used, then dustbin should be > 0
896
897	if ((microphys.and..not.doubleq).or.
898	& (microphys.and..not.water)) then
899	print*,'if microphys is used, then doubleq,'
900	print*,'and water must be used!'
901	call abort_physic(modname,
902	& "microphys requires water and doubleq",1)
903	endif
904	if (microphys.and..not.scavenging) then
905	print*,''
906	print*,'----------------WARNING-----------------'
907	print*,'microphys is used without scavenging !!!'
908	print*,'----------------WARNING-----------------'
909	print*,''
910	endif
911
912	if ((scavenging.and..not.microphys).or.
913	& (scavenging.and.(dustbin.lt.1)))then
914	print*,'if scavenging is used, then microphys'
915	print*,'must be used!'
916	call abort_physic(modname,
917	& "scavenging requires microphys",1)
918	endif
919
920	! Instantaneous scavenging by CO2
921	! -> expected to be replaced by scavenging with microphysics (flag scavenging) one day
922	write(,)"Dust scavenging by instantaneous CO2 snowfall ?"
923	scavco2cond=.false. ! default value
924	call getin_p("scavco2cond",scavco2cond)
925	write(,)" scavco2cond = ",scavco2cond
926	! Test of incompatibility:
927	! if scavco2cond is used, then dustbin should be > 0
928	if (scavco2cond.and.(dustbin.lt.1))then
929	print*,'if scavco2cond is used, then dustbin should be > 0'
930	call abort_physic(modname,
931	& "scavco2cond requires dustbin > 0",1)
932	endif
933	! if co2clouds is used, then there is no need for scavco2cond
934	if (co2clouds.and.scavco2cond) then
935	print*,''
936	print*,'----------------WARNING-----------------'
937	print*,' microphys scavenging is used so '
938	print*,' no need for scavco2cond !!! '
939	print*,'----------------WARNING-----------------'
940	print*,''
941	call abort_physic(modname,
942	& "incompatible co2cloud and scavco2cond options",1)
943	endif
944
945	! Test of incompatibility:
946
947	write(,) "Permanent water caps at poles ?",
948	& " .true. is RECOMMENDED"
949	write(,) "(with .true., North cap is a source of water ",
950	& "and South pole is a cold trap)"
951	caps=.true. ! default value
952	call getin_p("caps",caps)
953	write(,) " caps = ",caps
954
955	! Water ice albedo
956	write(,) "water cap albedo ?"
957	albedo_h2o_cap=0.35
958	call getin_p("albedo_h2o_cap",albedo_h2o_cap)
959	write(,) " albedo_h2o_cap = ",albedo_h2o_cap
960	write(,) "water frist albedo ?"
961	albedo_h2o_frost=0.35
962	call getin_p("albedo_h2o_frost",albedo_h2o_frost)
963	write(,) " albedo_h2o_frost = ",albedo_h2o_frost
964	write(,) "frost thickness threshold for albedo ?"
965	frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2)
966	call getin_p("frost_albedo_threshold",
967	& frost_albedo_threshold)
968	write(,) " frost_albedo_threshold = ",
969	& frost_albedo_threshold
970	write(,)"Water cap albedo is unchanged by water frost",
971	& " deposition (default is false)"
972	cst_cap_albedo=.false. ! default value
973	call getin_p("cst_cap_albedo",cst_cap_albedo)
974	write(,)"cst_cap_albedo = ",cst_cap_albedo
975
976	! Watercap evolution & refill (with discriminated albedos) (JN 2021)
977	write(,)"Watercap is replenished by water frost",
978	& " accumulation (default is false)"
979	refill_watercap=.false. ! default value
980	call getin_p("refill_watercap",refill_watercap)
981	write(,)"refill_watercap= ",refill_watercap
982	! frost thickness threshold for refill_watercap (ice metamorphism)
983	write(,) "frost thickness threshold for metamorphism ?",
984	& "ie converted into watercap",
985	& "only if refill_watercap is .true."
986	frost_metam_threshold=0.05 ! (i.e 0.05 kg.m-2)
987	call getin_p("frost_metam_threshold",
988	& frost_metam_threshold)
989	write(,) " frost_metam_threshold = ",
990	& frost_metam_threshold
991
992	! inert_h2o_ice
993	write(,) "water ice thermal inertia ?"
994	inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
995	call getin_p("inert_h2o_ice",inert_h2o_ice)
996	write(,) " inert_h2o_ice = ",inert_h2o_ice
997
998	! TMP: old_wsublimation_scheme
999	write(,) "Old water sublimation scheme?"
1000	old_wsublimation_scheme = .true.
1001	call getin_p("old_wsublimation_scheme",old_wsublimation_scheme)
1002	write(,) "old_wsublimation_scheme",old_wsublimation_scheme
1003
1004	! call Titus crocus line -- DEFAULT IS NONE
1005	write(,) "Titus crocus line ?"
1006	tituscap=.false. ! default value
1007	call getin_p("tituscap",tituscap)
1008	write(,) "tituscap",tituscap
1009
1010	! Chemistry:
1011	write(,) "photochemistry: include chemical species"
1012	photochem=.false. ! default value
1013	call getin_p("photochem",photochem)
1014	write(,) " photochem = ",photochem
1015
1016	write(,) "Compute chemistry (if photochem is .true.)",
1017	& "every ichemistry physics step (default: ichemistry=1)"
1018	ichemistry=1
1019	call getin_p("ichemistry",ichemistry)
1020	write(,) " ichemistry = ",ichemistry
1021
1022
1023	! SCATTERERS
1024	write(,) "how many scatterers?"
1025	naerkind=1 ! default value
1026	call getin_p("naerkind",naerkind)
1027	write(,)" naerkind = ",naerkind
1028
1029	! Test of incompatibility
1030	c Logical tests for radiatively active water-ice clouds:
1031	IF ( (activice.AND.(.NOT.water)).OR.
1032	& (activice.AND.(naerkind.LT.2)) ) THEN
1033	WRITE(,) 'If activice is TRUE, water has to be set'
1034	WRITE(,) 'to TRUE, and "naerkind" must be at least'
1035	WRITE(,) 'equal to 2.'
1036	call abort_physic(modname,
1037	& "radiatively active dust and water"//
1038	& " require naerkind > 1",1)
1039	ENDIF
1040
1041	!------------------------------------------
1042	!------------------------------------------
1043	! once naerkind is known allocate arrays
1044	! -- we do it here and not in phys_var_init
1045	! -- because we need to know naerkind first
1046	CALL ini_scatterers(ngrid,nlayer)
1047	!------------------------------------------
1048	!------------------------------------------
1049
1050
1051	c Please name the different scatterers here ----------------
1052	name_iaer(1) = "dust_conrath" !! default choice is good old Conrath profile
1053	IF (doubleq.AND.active) name_iaer(1) = "dust_doubleq" !! two-moment scheme
1054
1055	if (nq.gt.1) then
1056	! trick to avoid problems compiling with 1 tracer
1057	! and picky compilers who know name_iaer(2) is out of bounds
1058	j=2
1059	IF (rdstorm.AND..NOT.activice.AND..NOT.topflows) then
1060	name_iaer(j) = "stormdust_doubleq" !! storm dust two-moment scheme
1061	j = j+1
1062	END IF
1063
1064	IF (rdstorm.AND.water.AND.activice.AND..NOT.topflows) then
1065	name_iaer(j) = "stormdust_doubleq"
1066	j = j+1
1067	END IF
1068
1069	IF (topflows.AND..NOT.activice.AND..NOT.rdstorm) then
1070	name_iaer(j) = "topdust_doubleq" !! storm dust two-moment scheme
1071	j = j+1
1072	END IF
1073
1074	IF (topflows.AND.water.AND.activice.AND..NOT.rdstorm) then
1075	name_iaer(j) = "topdust_doubleq"
1076	j = j+1
1077	END IF
1078
1079	IF (rdstorm.AND.topflows.AND..NOT.activice) THEN
1080	name_iaer(j) = "stormdust_doubleq"
1081	name_iaer(j+1) = "topdust_doubleq"
1082	j = j+2
1083	ENDIF
1084
1085	IF (rdstorm.AND.topflows.AND.water.AND.activice) THEN
1086	name_iaer(j) = "stormdust_doubleq"
1087	name_iaer(j+1) = "topdust_doubleq"
1088	j = j+2
1089	ENDIF
1090
1091	IF (water.AND.activice) then
1092	name_iaer(j) = "h2o_ice" !! radiatively-active clouds
1093	j = j+1
1094	END IF
1095
1096	IF (co2clouds.AND.activeco2ice) then
1097	name_iaer(j) = "co2_ice" !! radiatively-active co2 clouds
1098	j = j+1
1099	ENDIF
1100
1101	IF (submicron.AND.active) then
1102	name_iaer(j) = "dust_submicron" !! JBM experimental stuff
1103	j = j+1
1104	ENDIF
1105	endif ! of if (nq.gt.1)
1106	c ----------------------------------------------------------
1107
1108	! Adsorption
1109	adsorption_soil = .false. ! default value
1110	call getin_p("adsorption_soil",adsorption_soil)
1111	if (adsorption_soil .and. (.not. water)) then
1112	write(,)"Adsorption can be run only if water = True"
1113	call abort_physic(modname,
1114	& "Adsorption must be used with water = true",1)
1115	endif
1116
1117	choice_ads = 0 ! default value (no adsorption)
1118	if (adsorption_soil) choice_ads = 1
1119	call getin_p("choice_ads",choice_ads)
1120
1121	ads_const_D = .false.
1122	call getin_p("ads_const_D",ads_const_D)
1123
1124	ads_massive_ice = .false.
1125	call getin_p("ads_massive_ice",ads_massive_ice)
1126
1127	poreice_tifeedback = .false.
1128	call getin_p("poreice_tifeedback",poreice_tifeedback)
1129	write(,) 'poreice_tifeedback=',poreice_tifeedback
1130	if (poreice_tifeedback .and. (.not. adsorption_soil)) then
1131	write(,)"Pore ice TI feedback can be run only if
1132	& adsorption_soil = True"
1133	call abort_physic(modname,
1134	& "Pore ice TI feedback be used with adsorption_soil
1135	& = true",1)
1136	endif
1137	c ----------------------------------------------------------
1138
1139	! THERMOSPHERE
1140
1141	write(,) "call thermosphere ?"
1142	callthermos=.false. ! default value
1143	call getin_p("callthermos",callthermos)
1144	write(,) " callthermos = ",callthermos
1145
1146
1147	write(,) " water included without cycle ",
1148	& "(only if water=.false.)"
1149	thermoswater=.false. ! default value
1150	call getin_p("thermoswater",thermoswater)
1151	write(,) " thermoswater = ",thermoswater
1152
1153	write(,) "call thermal conduction ?",
1154	& " (only if callthermos=.true.)"
1155	callconduct=.false. ! default value
1156	call getin_p("callconduct",callconduct)
1157	write(,) " callconduct = ",callconduct
1158
1159	write(,) "call EUV heating ?",
1160	& " (only if callthermos=.true.)"
1161	calleuv=.false. ! default value
1162	call getin_p("calleuv",calleuv)
1163	write(,) " calleuv = ",calleuv
1164
1165	write(,) "call molecular viscosity ?",
1166	& " (only if callthermos=.true.)"
1167	callmolvis=.false. ! default value
1168	call getin_p("callmolvis",callmolvis)
1169	write(,) " callmolvis = ",callmolvis
1170
1171	write(,) "call molecular diffusion ?",
1172	& " (only if callthermos=.true.)"
1173	callmoldiff=.false. ! default value
1174	call getin_p("callmoldiff",callmoldiff)
1175	write(,) " callmoldiff = ",callmoldiff
1176
1177
1178	write(,) "call thermospheric photochemistry ?",
1179	& " (only if callthermos=.true.)"
1180	thermochem=.false. ! default value
1181	call getin_p("thermochem",thermochem)
1182	write(,) " thermochem = ",thermochem
1183
1184	write(,) "Method to include solar variability"
1185	write(,) "0-> fixed value of E10.7 (fixed_euv_value); ",
1186	& "1-> daily evolution of E10.7 (for given solvaryear)"
1187	solvarmod=1
1188	call getin_p("solvarmod",solvarmod)
1189	write(,) " solvarmod = ",solvarmod
1190
1191	write(,) "Fixed euv (for solvarmod==0) 10.7 value?"
1192	write(,) " (min=80 , ave=140, max=320)"
1193	fixed_euv_value=140 ! default value
1194	call getin_p("fixed_euv_value",fixed_euv_value)
1195	write(,) " fixed_euv_value = ",fixed_euv_value
1196
1197	write(,) "Solar variability as observed for MY: "
1198	write(,) "Only if solvarmod=1"
1199	solvaryear=24
1200	call getin_p("solvaryear",solvaryear)
1201	write(,) " solvaryear = ",solvaryear
1202
1203	write(,) "UV heating efficiency:",
1204	& "measured values between 0.19 and 0.23 (Fox et al. 1996)",
1205	& "lower values may be used to compensate low 15 um cooling"
1206	euveff=0.21 !default value
1207	call getin_p("euveff",euveff)
1208	write(,) " euveff = ", euveff
1209
1210
1211	if (.not.callthermos) then
1212	if (thermoswater) then
1213	print*,'if thermoswater is set, callthermos must be true'
1214	call abort_physic(modname,
1215	& "thermoswater requires callthermos",1)
1216	endif
1217	if (callconduct) then
1218	print*,'if callconduct is set, callthermos must be true'
1219	call abort_physic(modname,
1220	& "callconduct requires callthermos",1)
1221	endif
1222	if (calleuv) then
1223	print*,'if calleuv is set, callthermos must be true'
1224	call abort_physic(modname,
1225	& "calleuv requires callthermos",1)
1226	endif
1227	if (callmolvis) then
1228	print*,'if callmolvis is set, callthermos must be true'
1229	call abort_physic(modname,
1230	& "callmolvis requires callthermos",1)
1231	endif
1232	if (callmoldiff) then
1233	print*,'if callmoldiff is set, callthermos must be true'
1234	call abort_physic(modname,
1235	& "callmoldiff requires callthermos",1)
1236	endif
1237	if (thermochem) then
1238	print*,'if thermochem is set, callthermos must be true'
1239	call abort_physic(modname,
1240	& "thermochem requires callthermos",1)
1241	endif
1242	endif
1243
1244	! Test of incompatibility:
1245	! if photochem is used, then water should be used too
1246
1247	if (photochem.and..not.water) then
1248	print*,'if photochem is used, water should be used too'
1249	call abort_physic(modname,
1250	& "photochem requires water",1)
1251	endif
1252
1253	! if callthermos is used, then thermoswater should be used too
1254	! (if water not used already)
1255
1256	if (callthermos .and. .not.water) then
1257	if (callthermos .and. .not.thermoswater) then
1258	print*,'if callthermos is used, water or thermoswater
1259	& should be used too'
1260	call abort_physic(modname,
1261	& "callthermos requires water or thermoswater",1)
1262	endif
1263	endif
1264
1265	PRINT*,'--------------------------------------------'
1266	PRINT*
1267	PRINT*
1268	ELSE
1269	write(,)
1270	write(,) 'Cannot read file callphys.def. Is it here ?'
1271	call abort_physic(modname,
1272	& "missing callphys.def file",1)
1273	ENDIF
1274
1275	8000 FORMAT(t5,a12,l8)
1276	8001 FORMAT(t5,a12,i8)
1277
1278	PRINT*
1279	PRINT*,'conf_phys: daysec',daysec
1280	PRINT*
1281	PRINT*,'conf_phys: The radiative transfer is computed:'
1282	PRINT*,' each ',iradia,' physical time-step'
1283	PRINT,' or each ',iradiadtphys,' seconds'
1284	PRINT*
1285	! --------------------------------------------------------------
1286	! Managing the Longwave radiative transfer
1287	! --------------------------------------------------------------
1288
1289	! In most cases, the run just use the following values :
1290	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1291	callemis=.true.
1292	! ilwd=10*int(daysec/dtphys) ! bug before 22/10/01
1293	ilwd=1
1294	ilwn=1 !2
1295	ilwb=1 !2
1296	linear=.true.
1297	ncouche=3
1298	alphan=0.4
1299	semi=0
1300
1301	! BUT people working hard on the LW may want to read them in 'radia.def'
1302	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1303	!$OMP MASTER
1304	OPEN(99,file='radia.def',status='old',form='formatted'
1305	. ,iostat=ierr)
1306	IF(ierr.EQ.0) THEN
1307	write(,) 'conf_phys: Reading radia.def !!!'
1308	READ(99,fmt='(a)') ch1
1309	READ(99,*) callemis
1310	WRITE(*,8000) ch1,callemis
1311
1312	READ(99,fmt='(a)') ch1
1313	READ(99,*) iradia
1314	WRITE(*,8001) ch1,iradia
1315
1316	READ(99,fmt='(a)') ch1
1317	READ(99,*) ilwd
1318	WRITE(*,8001) ch1,ilwd
1319
1320	READ(99,fmt='(a)') ch1
1321	READ(99,*) ilwn
1322	WRITE(*,8001) ch1,ilwn
1323
1324	READ(99,fmt='(a)') ch1
1325	READ(99,*) linear
1326	WRITE(*,8000) ch1,linear
1327
1328	READ(99,fmt='(a)') ch1
1329	READ(99,*) ncouche
1330	WRITE(*,8001) ch1,ncouche
1331
1332	READ(99,fmt='(a)') ch1
1333	READ(99,*) alphan
1334	WRITE(,) ch1,alphan
1335
1336	READ(99,fmt='(a)') ch1
1337	READ(99,*) ilwb
1338	WRITE(*,8001) ch1,ilwb
1339
1340
1341	READ(99,fmt='(a)') ch1
1342	READ(99,'(l1)') callg2d
1343	WRITE(*,8000) ch1,callg2d
1344
1345	READ(99,fmt='(a)') ch1
1346	READ(99,*) semi
1347	WRITE(,) ch1,semi
1348	end if
1349	CLOSE(99)
1350	!$OMP END MASTER
1351	#ifndef MESOSCALE
1352	call bcast(ch1)
1353	call bcast(callemis)
1354	call bcast(iradia)
1355	call bcast(ilwd)
1356	call bcast(ilwn)
1357	call bcast(linear)
1358	call bcast(ncouche)
1359	call bcast(alphan)
1360	call bcast(ilwb)
1361	call bcast(callg2d)
1362	call bcast(semi)
1363	#endif
1364
1365	END SUBROUTINE conf_phys
1366
1367	END MODULE conf_phys_mod

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

Download in other formats:

Original Format