Context Navigation

conf_phys.F @ 3468

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

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

Download in other formats:

Original Format