Context Navigation

conf_phys.F @ 2437

Last change on this file since 2437 was 2417, checked in by emillour, 5 years ago

Mars GCM:
Add a new scheme to handle the dust and its radiative impact. Triggered using
a new flag dustscaling_mode=2 (dustscaling_mod=0: no rescaling at all, and
dustscaling_mode=1: full rescaling using tauscaling, GCMv5.3 style). Rescaling
is then only done on the radiative impact (see dust_scaling_mod.F90) of dust.
Moreover the scaling factor "dust_rad_adjust" is evaluated using the target dust
scenario opacity for the next sol and left to evolve linearly until then to not
impose the diurnal evolution of dust.
In practice, main changes or additions in the code are:

renamed flag "tauscaling_mode" as "dustscaling_mode"
moved parameter "t_scenario_sol" to "dust_param_mod"
adapted "compute_dustscaling" routine in "dust_scaling_mod"
added module "dust_rad_adjust_mod"
2D fields "dust_rad_adjust_prev" and "dust_rad_adjust_next" required to compute coefficient "dust_rad_adjust" need to be stored in (re)startfi files

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

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

Download in other formats:

Original Format