Context Navigation

conf_phys.F @ 3807

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

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

Download in other formats:

Original Format