Context Navigation

conf_phys.F @ 3700

Last change on this file since 3700 was 3582, checked in by jbclement, 10 months ago

Mars PCM:

Correction of r3581: the key-words to get 'albedo_perennialco2(1:2)" from the "callphys.def" are now 'albedo_perennialco2_north' and 'albedo_perennialco2_south'.
Moving 'albedo_perennialco2' from the 'paleoclimate' module to the 'surfdat_h' module.

JBC

File size: 50.3 KB

Line
1	MODULE conf_phys_mod
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	SUBROUTINE conf_phys(ngrid,nlayer,nq)
8
9	!=======================================================================
10	!
11	! purpose:
12	! -------
13	!
14	! Initialisation for the physical parametrisations
15	! flags (i.e. run-time options) of the Mars PCM
16	!-----------------------------------------------------------------------
17
18	USE ioipsl_getin_p_mod, ONLY : getin_p
19	use tracer_mod, only : nuice_sed, ccn_factor, nuiceco2_sed,
20	& nuice_ref,nuiceco2_ref
21	use surfdat_h, only: albedo_h2o_cap,albedo_h2o_frost,
22	& albedo_perennialco2,
23	& frost_albedo_threshold, inert_h2o_ice,
24	& frost_metam_threshold,old_wsublimation_scheme
25	use time_phylmdz_mod, only: steps_per_sol,outputs_per_sol,iphysiq,
26	& ecritstart,daysec,dtphys
27	use dimradmars_mod, only: naerkind, name_iaer,
28	& ini_scatterers,tauvis
29	use datafile_mod, only: datadir
30	use wstats_mod, only: callstats
31	use writediagsoil_mod, only: diagsoil
32	use calchim_mod, only: ichemistry
33	use co2condens_mod, only: scavco2cond
34	use dust_param_mod, only: dustbin, doubleq, submicron, active,
35	& lifting, freedust, callddevil,
36	& dustscaling_mode,
37	& reff_driven_IRtoVIS_scenario
38	use aeropacity_mod, only: iddist, topdustref
39	USE mod_phys_lmdz_transfert_para, ONLY: bcast
40	USE paleoclimate_mod,ONLY: paleoclimate,
41	& lag_layer, include_waterbuoyancy
42	use microphys_h, only: mteta
43	use comsoil_h, only: adsorption_soil, choice_ads,ads_const_D,
44	& ads_massive_ice
45	use nonoro_gwd_mix_mod, only: calljliu_gwimix
46
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	real :: albedo_perennialco2_north, albedo_perennialco2_south
60
61	CHARACTER ch1*12
62	#ifndef MESOSCALE
63	! read in some parameters from "run.def" for physics,
64	! or shared between dynamics and physics.
65	ecritphy=-666 ! dummy default value
66	call getin_p("ecritphy",ecritphy) ! frequency of outputs in physics,
67	! in dynamical steps
68	if (ecritphy/=-666) then
69	call abort_physic(modname,
70	& "Error: parameter ecritphy is obsolete! Remove it! "//
71	& "And use outputs_per_sol instead",1)
72	endif
73
74	outputs_per_sol=4 ! default value, 4 outputs per sol
75	call getin_p("outputs_per_sol",outputs_per_sol)
76	! check that indeed it is possible for given physics time step
77	if (modulo(steps_per_sol,outputs_per_sol)/=0) then
78	write(,) "Error: outputs_per_sol = ",outputs_per_sol
79	write(,) " is not a divisor of number of steps per sol = ",
80	& steps_per_sol
81	call abort_physic(modname,
82	& "Error: inadequate value for outputs_per_sol",1)
83	endif
84
85	iphysiq=20 ! default value
86	call getin_p("iphysiq",iphysiq) ! call physics every iphysiq dyn step
87	ecritstart=0 ! default value
88	call getin_p("ecritstart",ecritstart) ! write a restart every ecristart steps
89	#endif
90
91	! --------------------------------------------------------------
92	! Reading the "callphys.def" file controlling some key options
93	! --------------------------------------------------------------
94
95	!$OMP MASTER
96	! check that 'callphys.def' file is around
97	OPEN(99,file='callphys.def',status='old',form='formatted'
98	& ,iostat=ierr)
99	CLOSE(99)
100	!$OMP END MASTER
101	call bcast(ierr)
102	! ierr=0
103
104	IF(ierr.EQ.0) THEN
105	PRINT*
106	PRINT*
107	PRINT*,'--------------------------------------------'
108	PRINT*,' conf_phys: Parameters for the physics (callphys.def)'
109	PRINT*,'--------------------------------------------'
110
111	write(,) "Directory where external input files are:"
112	! default path is set in datafile_mod
113	call getin_p("datadir",datadir)
114	write(,) " datadir = ",trim(datadir)
115
116	write(,) "Initialize physics with startfi.nc file ?"
117	startphy_file=.true.
118	call getin_p("startphy_file",startphy_file)
119	write(,) "startphy_file", startphy_file
120
121	write(,) "Diurnal cycle ?"
122	write(,) "(if diurnal=False, diurnal averaged solar heating)"
123	diurnal=.true. ! default value
124	call getin_p("diurnal",diurnal)
125	write(,) " diurnal = ",diurnal
126
127	write(,) "Seasonal cycle ?"
128	write(,) "(if season=False, Ls stays constant, to value ",
129	& "set in 'start'"
130	season=.true. ! default value
131	call getin_p("season",season)
132	write(,) " season = ",season
133
134	write(,) "Save statistics in file stats.nc ?"
135	#ifdef MESOSCALE
136	callstats=.false. ! default value
137	#else
138	callstats=.true. ! default value
139	#endif
140	call getin_p("callstats",callstats)
141	write(,) " callstats = ",callstats
142
143	write(,) "Write sub-surface fields in file diagsoil.nc ?"
144	diagsoil=.false. ! default value
145	call getin_p("diagsoil",diagsoil)
146	write(,) " diagsoil = ",diagsoil
147
148	write(,) "Save EOF profiles in file 'profiles' for ",
149	& "Climate Database?"
150	calleofdump=.false. ! default value
151	call getin_p("calleofdump",calleofdump)
152	write(,) " calleofdump = ",calleofdump
153
154	write(,) "Dust scenario: 1=constant dust (read from startfi",
155	& " or set as tauvis); 2=Viking scenario; =3 MGS scenario,",
156	& "=6 cold (low dust) scenario; =7 warm (high dust) scenario ",
157	& "=24,25 ... 30 :Mars Year 24, ... or 30 from TES assimilation"
158	iaervar=3 ! default value
159	call getin_p("iaervar",iaervar)
160	write(,) " iaervar = ",iaervar
161
162	write(,) "Reference (visible) dust opacity at 610 Pa ",
163	& "(matters only if iaervar=1)"
164	! NB: default value of tauvis is set/read in startfi.nc file
165	call getin_p("tauvis",tauvis)
166	write(,) " tauvis = ",tauvis
167
168	write(,) "Dust vertical distribution:"
169	write(,) "(=1 top set by topdustref parameter;",
170	& " =2 Viking scenario; =3 MGS scenario)"
171	iddist=3 ! default value
172	call getin_p("iddist",iddist)
173	write(,) " iddist = ",iddist
174
175	write(,) "Dust top altitude (km). (Matters only if iddist=1)"
176	topdustref= 90.0 ! default value
177	call getin_p("topdustref",topdustref)
178	write(,) " topdustref = ",topdustref
179
180	write(,) "Prescribed surface thermal flux (H/(rho*cp),K m/s)"
181	tke_heat_flux=0. ! default value
182	call getin_p("tke_heat_flux",tke_heat_flux)
183	write(,) " tke_heat_flux = ",tke_heat_flux
184	write(,) " 0 means the usual schemes are computing"
185
186	write(,) "call radiative transfer ?"
187	callrad=.true. ! default value
188	call getin_p("callrad",callrad)
189	write(,) " callrad = ",callrad
190
191	write(,) "call slope insolation scheme ?",
192	& "(matters only if callrad=T)"
193	#ifdef MESOSCALE
194	callslope=.true. ! default value
195	#else
196	callslope=.false. ! default value (not supported yet)
197	#endif
198	call getin_p("callslope",callslope)
199	write(,) " callslope = ",callslope
200
201	write(,) "call NLTE radiative schemes ?",
202	& "(matters only if callrad=T)"
203	callnlte=.false. ! default value
204	call getin_p("callnlte",callnlte)
205	write(,) " callnlte = ",callnlte
206
207	nltemodel=0 !default value
208	write(,) "NLTE model?"
209	write(,) "0 -> old model, static O"
210	write(,) "1 -> old model, dynamic O"
211	write(,) "2 -> new model"
212	write(,) "(matters only if callnlte=T)"
213	call getin_p("nltemodel",nltemodel)
214	write(,) " nltemodel = ",nltemodel
215
216	write(,) "call CO2 NIR absorption ?",
217	& "(matters only if callrad=T)"
218	callnirco2=.false. ! default value
219	call getin_p("callnirco2",callnirco2)
220	write(,) " callnirco2 = ",callnirco2
221
222	write(,) "New NIR NLTE correction ?",
223	$ "0-> old model (no correction)",
224	$ "1-> new correction",
225	$ "(matters only if callnirco2=T)"
226	#ifdef MESOSCALE
227	nircorr=0 !default value. this is OK below 60 km.
228	#else
229	nircorr=0 !default value
230	#endif
231	call getin_p("nircorr",nircorr)
232	write(,) " nircorr = ",nircorr
233
234	write(,) "call turbulent vertical diffusion ?"
235	calldifv=.true. ! default value
236	call getin_p("calldifv",calldifv)
237	write(,) " calldifv = ",calldifv
238
239	write(,) "call thermals ?"
240	calltherm=.false. ! default value
241	call getin_p("calltherm",calltherm)
242	write(,) " calltherm = ",calltherm
243
244	write(,) "call convective adjustment ?"
245	calladj=.true. ! default value
246	call getin_p("calladj",calladj)
247	write(,) " calladj = ",calladj
248
249	if (calltherm .and. calladj) then
250	print*,'!!! PLEASE NOTE !!!'
251	print*,'convective adjustment is on'
252	print*,'but since thermal plume model is on'
253	print*,'convadj is only activated above the PBL'
254	endif
255
256	write(,) "used latest version of yamada scheme?"
257	callyamada4=.true. ! default value
258	call getin_p("callyamada4",callyamada4)
259	write(,) " callyamada4 = ",callyamada4
260
261	write(,) "used latest version of ATKE scheme?"
262	callatke =.false. ! default value
263	call getin_p("callatke",callatke)
264	write(,) " callatke = ",callatke
265
266	if(callyamada4.and.callatke) then
267	print*,"You can't use both yamada and atke"
268	print*,"Choose either Yamada or ATKE"
269	call abort_physic(modname,
270	& "Can't use both yamada and ATKE",1)
271	endif
272
273	if (calltherm .and. .not.callyamada4) then
274	print*,'!!!! WARNING WARNING WARNING !!!!'
275	print*,'if calltherm=T we strongly advise that '
276	print*,'you set the flag callyamada4 to T '
277	print*,'!!!! WARNING WARNING WARNING !!!!'
278	endif
279
280	write(,) "call Richardson-based surface layer ?"
281	callrichsl=.false. ! default value
282	call getin_p("callrichsl",callrichsl)
283	write(,) " callrichsl = ",callrichsl
284
285	if (calltherm .and. .not.callrichsl) then
286	print*,'WARNING WARNING WARNING'
287	print*,'if calltherm=T we strongly advise that '
288	print*,'you use the new surface layer scheme '
289	print*,'by setting callrichsl=T '
290	endif
291
292	if (calladj .and. callrichsl .and. (.not. calltherm)) then
293	print*,'You should not be calling the convective adjustment
294	& scheme with the Richardson surface-layer and without the thermals
295	&. This approach is not
296	& physically consistent and can lead to unrealistic friction
297	& values.'
298	print*,'If you want to use the Ri. surface-layer, either
299	& activate thermals OR de-activate the convective adjustment.'
300	call abort_physic(modname,
301	& "Richardson layer must be used with thermals",1)
302	endif
303
304	write(,) "call CO2 condensation ?"
305	callcond=.true. ! default value
306	call getin_p("callcond",callcond)
307	write(,) " callcond = ",callcond
308
309	write(,)"call thermal conduction in the soil ?"
310	callsoil=.true. ! default value
311	call getin_p("callsoil",callsoil)
312	write(,) " callsoil = ",callsoil
313
314
315	write(,)"call Lott's gravity wave/subgrid topography ",
316	& "scheme ?"
317	calllott=.true. ! default value
318	call getin_p("calllott",calllott)
319	write(,)" calllott = ",calllott
320
321	write(,)"call Lott's non-oro GWs parameterisation ",
322	& "scheme ?"
323	calllott_nonoro=.false. ! default value
324	call getin_p("calllott_nonoro",calllott_nonoro)
325	write(,)" calllott_nonoro = ",calllott_nonoro
326
327	write(,)"call jliu's nonoro GWs-induced mixing ",
328	& "scheme ?"
329	calljliu_gwimix=.false. ! default value
330	call getin_p("calljliu_gwimix",calljliu_gwimix)
331	write(,)" calljliu_gwimix = ",calljliu_gwimix
332
333
334	! rocket dust storm injection scheme
335	write(,)"call rocket dust storm parametrization"
336	rdstorm=.false. ! default value
337	call getin_p("rdstorm",rdstorm)
338	write(,)" rdstorm = ",rdstorm
339	! rocket dust storm detrainment coefficient
340	coeff_detrainment=0.02 ! default value
341	call getin_p("coeff_detrainment",coeff_detrainment)
342	write(,)" coeff_detrainment = ",coeff_detrainment
343
344	! entrainment by mountain top dust flows scheme
345	write(,)"call mountain top dust flows parametrization"
346	topflows=.false. ! default value
347	call getin_p("topflows",topflows)
348	write(,)" topflows = ",topflows
349
350	! latent heat release from ground water ice sublimation/condensation
351	write(,)"latent heat release during sublimation",
352	& " /condensation of ground water ice"
353	latentheat_surfwater=.true. ! default value
354	call getin_p("latentheat_surfwater",latentheat_surfwater)
355	write(,)" latentheat_surfwater = ",latentheat_surfwater
356
357	write(,)"rad.transfer is computed every iradia",
358	& " physical timestep"
359	iradia=1 ! default value
360	call getin_p("iradia",iradia)
361	write(,)" iradia = ",iradia
362
363
364	write(,)"Output of the exchange coefficient mattrix ?",
365	& "(for diagnostics only)"
366	callg2d=.false. ! default value
367	call getin_p("callg2d",callg2d)
368	write(,)" callg2d = ",callg2d
369
370	write(,)"Rayleigh scattering : (should be .false. for now)"
371	rayleigh=.false.
372	call getin_p("rayleigh",rayleigh)
373	write(,)" rayleigh = ",rayleigh
374
375	! PALEOCLIMATE
376
377	write(,)"Using lag layer??"
378	lag_layer=.false.
379	call getin_p("lag_layer",lag_layer)
380	write(,) "lag_layer = ", lag_layer
381
382	write(,)"Is it paleoclimate run?"
383	paleoclimate=.false. ! default value
384	call getin_p("paleoclimate",paleoclimate)
385	write(,)"paleoclimate = ",paleoclimate
386
387	write(,)"Albedo for perennial CO2 ice (northern hemisphere)?"
388	albedo_perennialco2_north = 0.6 ! default value
389	call getin_p("albedo_perennialco2_north",albedo_perennialco2_
390	& north)
391	albedo_perennialco2(1) = albedo_perennialco2_north
392	write(,)"albedo_perennialco2(1) = ",albedo_perennialco2(1)
393
394	write(,)"Albedo for perennial CO2 ice (southern hemisphere)?"
395	albedo_perennialco2_south = 0.85 ! default value
396	call getin_p("albedo_perennialco2_south",albedo_perennialco2_
397	& south)
398	albedo_perennialco2(2) = albedo_perennialco2_south
399	write(,)"albedo_perennialco2(2) = ",albedo_perennialco2(2)
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