Context Navigation

conf_phys.F @ 2325

Last change on this file since 2325 was 2312, checked in by lrossi, 5 years ago

MARS GCM:
Implementation of HDO in the physics
New tracers hdo_vap and hdo_ice are added. Cf. traceur.def.isotopes in deftank for exemple of traceur.def.
All HDO related computations are activated by flag hdo=.true. in callphys.def. (see callphys.def.hdo in deftank for an example)
Additional option hdofrac (true by default) allows for turning on/off fractionation (for tests).
For now, only functional with simplified cloud scheme (so microphys=.false. and activice=.false. also recommended).
Initialisation of start files can be done with option 'inihdo' in newstart.
LR

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

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

Download in other formats:

Original Format