Context Navigation

conf_phys.F @ 2616

Last change on this file since 2616 was 2561, checked in by jnaar, 3 years ago

Water cycle flag update & addition :

"cap_albedo" renamed "cst_cap_albedo" (default false) : if true, water ice cap albedo remains constant even when frost with higher albedo condensates on it
"refill_watercap" added (default false) : turns h2o_ice_s into watercap when above a given threshold
"frost_metam_threshold" added (default 0.05 km.m-2) : threshold used by "refill_watercap"

File size: 43.1 KB

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

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

Download in other formats:

Original Format