Context Navigation

conf_phys.F @ 2511

Last change on this file since 2511 was 2508, checked in by jnaar, 4 years ago
The water ice albedo is now conceptually decoupled between old perennial ice and seasonal frost. For now, "old ice" is only on watercaptag = .true. The variable albedo_h2o_ice has been decomposed into albedo_h2o_cap and albedo_h2o_frost. Default values are both 0.35 and retrocompatible with old callphys.def and albedo_h2o_ice if needed. JN
File size: 41.4 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	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	! Activeco2ice
568	write(,) "Radiatively active transported atmospheric ",
569	& "Co2 ice ?"
570	activeco2ice=.false. ! default value
571	call getin_p("activeco2ice",activeco2ice)
572	write(,) " activeco2ice = ",activeco2ice
573
574	! sub-grid cloud fraction: fixed
575	write(,) "Fixed cloud fraction?"
576	CLFfixval=1.0 ! default value
577	call getin_p("CLFfixval",CLFfixval)
578	write(,) "CLFfixval=",CLFfixval
579	! sub-grid cloud fraction: varying
580	write(,) "Use partial nebulosity?"
581	CLFvarying=.false. ! default value
582	call getin_p("CLFvarying",CLFvarying)
583	write(,)"CLFvarying=",CLFvarying
584
585	!CO2 clouds scheme?
586	write(,) "Compute CO2 clouds (implies microphysical scheme)?"
587	co2clouds=.false. ! default value
588	call getin_p("co2clouds",co2clouds)
589	write(,) " co2clouds = ",co2clouds
590	!Can water ice particles serve as CCN for CO2clouds
591	write(,) "Use water ice as CO2 clouds CCN ?"
592	co2useh2o=.false. ! default value
593	call getin_p("co2useh2o",co2useh2o)
594	write(,) " co2useh2o = ",co2useh2o
595	!Do we allow a supply of meteoritic paricles to serve as CO2 ice CCN?
596	write(,) "Supply meteoritic particle for CO2 clouds ?"
597	meteo_flux=.false. !Default value
598	call getin_p("meteo_flux",meteo_flux)
599	write(,) " meteo_flux = ",meteo_flux
600	!Do we allow a sub-grid temperature distribution for the CO2 microphysics
601	write(,) "sub-grid temperature distribution for CO2 clouds?"
602	CLFvaryingCO2=.false. !Default value
603	call getin_p("CLFvaryingCO2",CLFvaryingCO2)
604	write(,) " CLFvaryingCO2 = ",CLFvaryingCO2
605	!Amplitude of the sub-grid temperature distribution for the CO2 microphysics
606	write(,) "sub-grid temperature amplitude for CO2 clouds?"
607	spantCO2=0 !Default value
608	call getin_p("spantCO2",spantCO2)
609	write(,) " spantCO2 = ",spantCO2
610	!Do you want to filter the sub-grid T distribution by a Saturation index?
611	write(,) "filter sub-grid temperature by Saturation index?"
612	satindexco2=.true.
613	call getin_p("satindexco2",satindexco2)
614	write(,) " satindexco2 = ",satindexco2
615
616
617	! thermal inertia feedback
618	write(,) "Activate the thermal inertia feedback ?"
619	tifeedback=.false. ! default value
620	call getin_p("tifeedback",tifeedback)
621	write(,) " tifeedback = ",tifeedback
622
623	! Test of incompatibility:
624
625	if (tifeedback.and..not.water) then
626	print*,'if tifeedback is used,'
627	print*,'water should be used too'
628	call abort_physic(modname,
629	& "tifeedback requires water",1)
630	endif
631
632	if (tifeedback.and..not.callsoil) then
633	print*,'if tifeedback is used,'
634	print*,'callsoil should be used too'
635	call abort_physic(modname,
636	& "tifeedback requires callsoil",1)
637	endif
638
639	if (activice.and..not.water) then
640	print*,'if activice is used, water should be used too'
641	call abort_physic(modname,
642	& "activeice requires water",1)
643	endif
644
645	if (water.and..not.tracer) then
646	print*,'if water is used, tracer should be used too'
647	call abort_physic(modname,
648	& "water requires tracer",1)
649	endif
650
651	if (hdo.and..not.water) then
652	print*,'if hdo is used, water should be used too'
653	call abort_physic(modname,
654	& "hd2 requires tracer",1)
655	endif
656
657
658	if (activeco2ice.and..not.co2clouds) then
659	print*,'if activeco2ice is used, co2clouds should be used too'
660	call abort_physic(modname,
661	& "activeco2ice requires co2clouds",1)
662	endif
663
664	! water ice clouds effective variance distribution for sedimentaion
665	write(,) "Sed effective variance for water ice clouds ?"
666	nuice_sed=0.45
667	call getin_p("nuice_sed",nuice_sed)
668	write(,) "water_param nueff Sedimentation:", nuice_sed
669
670	write(,) "Sed effective variance for CO2 clouds ?"
671	nuiceco2_sed=0.45
672	call getin_p("nuiceco2_sed",nuiceco2_sed)
673	write(,) "CO2 nueff Sedimentation:", nuiceco2_sed
674
675	write(,) "REF effective variance for CO2 clouds ?"
676	nuiceco2_ref=0.45
677	call getin_p("nuiceco2_ref",nuiceco2_ref)
678	write(,) "CO2 nueff Sedimentation:", nuiceco2_ref
679
680	write(,) "REF effective variance for water clouds ?"
681	nuice_ref=0.45
682	call getin_p("nuice_ref",nuice_ref)
683	write(,) "CO2 nueff Sedimentation:", nuice_ref
684
685
686	! ccn factor if no scavenging
687	write(,) "water param CCN reduc. factor ?"
688	ccn_factor = 4.5
689	call getin_p("ccn_factor",ccn_factor)
690	write(,)" ccn_factor = ",ccn_factor
691	write(,)"Careful: only used when microphys=F, otherwise"
692	write(,)"the contact parameter is used instead;"
693
694	! microphys
695	write(,)"Microphysical scheme for water-ice clouds?"
696	microphys=.false. ! default value
697	call getin_p("microphys",microphys)
698	write(,)" microphys = ",microphys
699
700	! supersat
701	write(,)"Allow super-saturation of water vapor?"
702	supersat=.true. ! default value
703	call getin_p("supersat",supersat)
704	write(,)"supersat = ",supersat
705
706	! microphysical parameter contact
707	write(,) "water contact parameter ?"
708	mteta = 0.95
709	call getin_p("mteta",mteta)
710	write(,) "mteta = ", mteta
711
712	! scavenging
713	write(,)"Dust scavenging by H2O/CO2 snowfall ?"
714	scavenging=.false. ! default value
715	call getin_p("scavenging",scavenging)
716	write(,)" scavenging = ",scavenging
717
718
719	! Test of incompatibility:
720	! if scavenging is used, then dustbin should be > 0
721
722	if ((microphys.and..not.doubleq).or.
723	& (microphys.and..not.water)) then
724	print*,'if microphys is used, then doubleq,'
725	print*,'and water must be used!'
726	call abort_physic(modname,
727	& "microphys requires water and doubleq",1)
728	endif
729	if (microphys.and..not.scavenging) then
730	print*,''
731	print*,'----------------WARNING-----------------'
732	print*,'microphys is used without scavenging !!!'
733	print*,'----------------WARNING-----------------'
734	print*,''
735	endif
736
737	if ((scavenging.and..not.microphys).or.
738	& (scavenging.and.(dustbin.lt.1)))then
739	print*,'if scavenging is used, then microphys'
740	print*,'must be used!'
741	call abort_physic(modname,
742	& "scavenging requires microphys",1)
743	endif
744
745	! Instantaneous scavenging by CO2
746	! -> expected to be replaced by scavenging with microphysics (flag scavenging) one day
747	write(,)"Dust scavenging by instantaneous CO2 snowfall ?"
748	scavco2cond=.false. ! default value
749	call getin_p("scavco2cond",scavco2cond)
750	write(,)" scavco2cond = ",scavco2cond
751	! Test of incompatibility:
752	! if scavco2cond is used, then dustbin should be > 0
753	if (scavco2cond.and.(dustbin.lt.1))then
754	print*,'if scavco2cond is used, then dustbin should be > 0'
755	call abort_physic(modname,
756	& "scavco2cond requires dustbin > 0",1)
757	endif
758	! if co2clouds is used, then there is no need for scavco2cond
759	if (co2clouds.and.scavco2cond) then
760	print*,''
761	print*,'----------------WARNING-----------------'
762	print*,' microphys scavenging is used so '
763	print*,' no need for scavco2cond !!! '
764	print*,'----------------WARNING-----------------'
765	print*,''
766	call abort_physic(modname,
767	& "incompatible co2cloud and scavco2cond options",1)
768	endif
769
770	! Test of incompatibility:
771
772	write(,) "Permanent water caps at poles ?",
773	& " .true. is RECOMMENDED"
774	write(,) "(with .true., North cap is a source of water ",
775	& "and South pole is a cold trap)"
776	caps=.true. ! default value
777	call getin_p("caps",caps)
778	write(,) " caps = ",caps
779
780	! JN : now separated between albedo_h2o_cap and
781	! albedo_h2o_frost. Retrocompatible with old
782	! callphys.def with albedo_h2o_ice
783	write(,) "water ice albedo ?"
784	albedo_h2o_cap=0.35
785	albedo_h2o_frost=0.35
786	call getin_p("albedo_h2o_ice",albedo_h2o_cap)
787	albedo_h2o_frost=albedo_h2o_cap
788	call getin_p("albedo_h2o_cap",albedo_h2o_cap)
789	write(,) " albedo_h2o_cap = ",albedo_h2o_cap
790	call getin_p("albedo_h2o_frost",albedo_h2o_frost)
791	write(,) " albedo_h2o_frost = ",albedo_h2o_frost
792	! inert_h2o_ice
793	write(,) "water ice thermal inertia ?"
794	inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
795	call getin_p("inert_h2o_ice",inert_h2o_ice)
796	write(,) " inert_h2o_ice = ",inert_h2o_ice
797	! frost_albedo_threshold
798	write(,) "frost thickness threshold for albedo ?"
799	frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2)
800	call getin_p("frost_albedo_threshold",
801	& frost_albedo_threshold)
802	write(,) " frost_albedo_threshold = ",
803	& frost_albedo_threshold
804
805	! call Titus crocus line -- DEFAULT IS NONE
806	write(,) "Titus crocus line ?"
807	tituscap=.false. ! default value
808	call getin_p("tituscap",tituscap)
809	write(,) "tituscap",tituscap
810
811	! Chemistry:
812	write(,) "photochemistry: include chemical species"
813	photochem=.false. ! default value
814	call getin_p("photochem",photochem)
815	write(,) " photochem = ",photochem
816
817	write(,) "Compute chemistry (if photochem is .true.)",
818	& "every ichemistry physics step (default: ichemistry=1)"
819	ichemistry=1
820	call getin_p("ichemistry",ichemistry)
821	write(,) " ichemistry = ",ichemistry
822
823
824	! SCATTERERS
825	write(,) "how many scatterers?"
826	naerkind=1 ! default value
827	call getin_p("naerkind",naerkind)
828	write(,)" naerkind = ",naerkind
829
830	! Test of incompatibility
831	c Logical tests for radiatively active water-ice clouds:
832	IF ( (activice.AND.(.NOT.water)).OR.
833	& (activice.AND.(naerkind.LT.2)) ) THEN
834	WRITE(,) 'If activice is TRUE, water has to be set'
835	WRITE(,) 'to TRUE, and "naerkind" must be at least'
836	WRITE(,) 'equal to 2.'
837	call abort_physic(modname,
838	& "radiatively active dust and water"//
839	& " require naerkind > 1",1)
840	ENDIF
841
842	!------------------------------------------
843	!------------------------------------------
844	! once naerkind is known allocate arrays
845	! -- we do it here and not in phys_var_init
846	! -- because we need to know naerkind first
847	CALL ini_scatterers(ngrid,nlayer)
848	!------------------------------------------
849	!------------------------------------------
850
851
852	c Please name the different scatterers here ----------------
853	name_iaer(1) = "dust_conrath" !! default choice is good old Conrath profile
854	IF (doubleq.AND.active) name_iaer(1) = "dust_doubleq" !! two-moment scheme
855
856	if (nq.gt.1) then
857	! trick to avoid problems compiling with 1 tracer
858	! and picky compilers who know name_iaer(2) is out of bounds
859	j=2
860	IF (rdstorm.AND..NOT.activice.AND..NOT.slpwind) then
861	name_iaer(j) = "stormdust_doubleq" !! storm dust two-moment scheme
862	j = j+1
863	END IF
864
865	IF (rdstorm.AND.water.AND.activice.AND..NOT.slpwind) then
866	name_iaer(j) = "stormdust_doubleq"
867	j = j+1
868	END IF
869
870	IF (slpwind.AND..NOT.activice.AND..NOT.rdstorm) then
871	name_iaer(j) = "topdust_doubleq" !! storm dust two-moment scheme
872	j = j+1
873	END IF
874
875	IF (slpwind.AND.water.AND.activice.AND..NOT.rdstorm) then
876	name_iaer(j) = "topdust_doubleq"
877	j = j+1
878	END IF
879
880	IF (rdstorm.AND.slpwind.AND..NOT.activice) THEN
881	name_iaer(j) = "stormdust_doubleq"
882	name_iaer(j+1) = "topdust_doubleq"
883	j = j+2
884	ENDIF
885
886	IF (rdstorm.AND.slpwind.AND.water.AND.activice) THEN
887	name_iaer(j) = "stormdust_doubleq"
888	name_iaer(j+1) = "topdust_doubleq"
889	j = j+2
890	ENDIF
891
892	IF (water.AND.activice) then
893	name_iaer(j) = "h2o_ice" !! radiatively-active clouds
894	j = j+1
895	END IF
896
897	IF (co2clouds.AND.activeco2ice) then
898	name_iaer(j) = "co2_ice" !! radiatively-active co2 clouds
899	j = j+1
900	ENDIF
901
902	IF (submicron.AND.active) then
903	name_iaer(j) = "dust_submicron" !! JBM experimental stuff
904	j = j+1
905	ENDIF
906	endif ! of if (nq.gt.1)
907	c ----------------------------------------------------------
908
909	! THERMOSPHERE
910
911	write(,) "call thermosphere ?"
912	callthermos=.false. ! default value
913	call getin_p("callthermos",callthermos)
914	write(,) " callthermos = ",callthermos
915
916
917	write(,) " water included without cycle ",
918	& "(only if water=.false.)"
919	thermoswater=.false. ! default value
920	call getin_p("thermoswater",thermoswater)
921	write(,) " thermoswater = ",thermoswater
922
923	write(,) "call thermal conduction ?",
924	& " (only if callthermos=.true.)"
925	callconduct=.false. ! default value
926	call getin_p("callconduct",callconduct)
927	write(,) " callconduct = ",callconduct
928
929	write(,) "call EUV heating ?",
930	& " (only if callthermos=.true.)"
931	calleuv=.false. ! default value
932	call getin_p("calleuv",calleuv)
933	write(,) " calleuv = ",calleuv
934
935	write(,) "call molecular viscosity ?",
936	& " (only if callthermos=.true.)"
937	callmolvis=.false. ! default value
938	call getin_p("callmolvis",callmolvis)
939	write(,) " callmolvis = ",callmolvis
940
941	write(,) "call molecular diffusion ?",
942	& " (only if callthermos=.true.)"
943	callmoldiff=.false. ! default value
944	call getin_p("callmoldiff",callmoldiff)
945	write(,) " callmoldiff = ",callmoldiff
946
947
948	write(,) "call thermospheric photochemistry ?",
949	& " (only if callthermos=.true.)"
950	thermochem=.false. ! default value
951	call getin_p("thermochem",thermochem)
952	write(,) " thermochem = ",thermochem
953
954	write(,) "Method to include solar variability"
955	write(,) "0-> fixed value of E10.7 (fixed_euv_value); ",
956	& "1-> daily evolution of E10.7 (for given solvaryear)"
957	solvarmod=1
958	call getin_p("solvarmod",solvarmod)
959	write(,) " solvarmod = ",solvarmod
960
961	write(,) "Fixed euv (for solvarmod==0) 10.7 value?"
962	write(,) " (min=80 , ave=140, max=320)"
963	fixed_euv_value=140 ! default value
964	call getin_p("fixed_euv_value",fixed_euv_value)
965	write(,) " fixed_euv_value = ",fixed_euv_value
966
967	write(,) "Solar variability as observed for MY: "
968	write(,) "Only if solvarmod=1"
969	solvaryear=24
970	call getin_p("solvaryear",solvaryear)
971	write(,) " solvaryear = ",solvaryear
972
973	write(,) "UV heating efficiency:",
974	& "measured values between 0.19 and 0.23 (Fox et al. 1996)",
975	& "lower values may be used to compensate low 15 um cooling"
976	euveff=0.21 !default value
977	call getin_p("euveff",euveff)
978	write(,) " euveff = ", euveff
979
980
981	if (.not.callthermos) then
982	if (thermoswater) then
983	print*,'if thermoswater is set, callthermos must be true'
984	call abort_physic(modname,
985	& "thermoswater requires callthermos",1)
986	endif
987	if (callconduct) then
988	print*,'if callconduct is set, callthermos must be true'
989	call abort_physic(modname,
990	& "callconduct requires callthermos",1)
991	endif
992	if (calleuv) then
993	print*,'if calleuv is set, callthermos must be true'
994	call abort_physic(modname,
995	& "calleuv requires callthermos",1)
996	endif
997	if (callmolvis) then
998	print*,'if callmolvis is set, callthermos must be true'
999	call abort_physic(modname,
1000	& "callmolvis requires callthermos",1)
1001	endif
1002	if (callmoldiff) then
1003	print*,'if callmoldiff is set, callthermos must be true'
1004	call abort_physic(modname,
1005	& "callmoldiff requires callthermos",1)
1006	endif
1007	if (thermochem) then
1008	print*,'if thermochem is set, callthermos must be true'
1009	call abort_physic(modname,
1010	& "thermochem requires callthermos",1)
1011	endif
1012	endif
1013
1014	! Test of incompatibility:
1015	! if photochem is used, then water should be used too
1016
1017	if (photochem.and..not.water) then
1018	print*,'if photochem is used, water should be used too'
1019	call abort_physic(modname,
1020	& "photochem requires water",1)
1021	endif
1022
1023	! if callthermos is used, then thermoswater should be used too
1024	! (if water not used already)
1025
1026	if (callthermos .and. .not.water) then
1027	if (callthermos .and. .not.thermoswater) then
1028	print*,'if callthermos is used, water or thermoswater
1029	& should be used too'
1030	call abort_physic(modname,
1031	& "callthermos requires water or thermoswater",1)
1032	endif
1033	endif
1034
1035	PRINT*,'--------------------------------------------'
1036	PRINT*
1037	PRINT*
1038	ELSE
1039	write(,)
1040	write(,) 'Cannot read file callphys.def. Is it here ?'
1041	call abort_physic(modname,
1042	& "missing callphys.def file",1)
1043	ENDIF
1044
1045	8000 FORMAT(t5,a12,l8)
1046	8001 FORMAT(t5,a12,i8)
1047
1048	PRINT*
1049	PRINT*,'conf_phys: daysec',daysec
1050	PRINT*
1051	PRINT*,'conf_phys: The radiative transfer is computed:'
1052	PRINT*,' each ',iradia,' physical time-step'
1053	PRINT,' or each ',iradiadtphys,' seconds'
1054	PRINT*
1055	! --------------------------------------------------------------
1056	! Managing the Longwave radiative transfer
1057	! --------------------------------------------------------------
1058
1059	! In most cases, the run just use the following values :
1060	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1061	callemis=.true.
1062	! ilwd=10*int(daysec/dtphys) ! bug before 22/10/01
1063	ilwd=1
1064	ilwn=1 !2
1065	ilwb=1 !2
1066	linear=.true.
1067	ncouche=3
1068	alphan=0.4
1069	semi=0
1070
1071	! BUT people working hard on the LW may want to read them in 'radia.def'
1072	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1073
1074	OPEN(99,file='radia.def',status='old',form='formatted'
1075	. ,iostat=ierr)
1076	IF(ierr.EQ.0) THEN
1077	write(,) 'conf_phys: Reading radia.def !!!'
1078	READ(99,fmt='(a)') ch1
1079	READ(99,*) callemis
1080	WRITE(*,8000) ch1,callemis
1081
1082	READ(99,fmt='(a)') ch1
1083	READ(99,*) iradia
1084	WRITE(*,8001) ch1,iradia
1085
1086	READ(99,fmt='(a)') ch1
1087	READ(99,*) ilwd
1088	WRITE(*,8001) ch1,ilwd
1089
1090	READ(99,fmt='(a)') ch1
1091	READ(99,*) ilwn
1092	WRITE(*,8001) ch1,ilwn
1093
1094	READ(99,fmt='(a)') ch1
1095	READ(99,*) linear
1096	WRITE(*,8000) ch1,linear
1097
1098	READ(99,fmt='(a)') ch1
1099	READ(99,*) ncouche
1100	WRITE(*,8001) ch1,ncouche
1101
1102	READ(99,fmt='(a)') ch1
1103	READ(99,*) alphan
1104	WRITE(,) ch1,alphan
1105
1106	READ(99,fmt='(a)') ch1
1107	READ(99,*) ilwb
1108	WRITE(*,8001) ch1,ilwb
1109
1110
1111	READ(99,fmt='(a)') ch1
1112	READ(99,'(l1)') callg2d
1113	WRITE(*,8000) ch1,callg2d
1114
1115	READ(99,fmt='(a)') ch1
1116	READ(99,*) semi
1117	WRITE(,) ch1,semi
1118	end if
1119	CLOSE(99)
1120
1121	END

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

Download in other formats:

Original Format