Context Navigation

conf_phys.F @ 3369

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

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

Download in other formats:

Original Format