Context Navigation

conf_phys.F @ 4008

Last change on this file since 4008 was 4008, checked in by aslmd, 4 weeks ago
MESOSCALE: use precompiling flags to hide instructions related to parallel computations (we consider physics as being like a 1D model without any attached dynamical core when we compile).
File size: 52.4 KB

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

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

Download in other formats:

Original Format