Context Navigation

conf_phys.F @ 2398

Last change on this file since 2398 was 2398, checked in by emillour, 5 years ago
Mars GCM: Some code cleanup: use "call abort_physic()" instead of "stop" or "call abort" EM
File size: 39.6 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	call abort_physic(modname,
635	& "hd2 requires tracer",1)
636	endif
637
638
639	! water ice clouds effective variance distribution for sedimentaion
640	write(,) "Sed effective variance for water ice clouds ?"
641	nuice_sed=0.45
642	call getin_p("nuice_sed",nuice_sed)
643	write(,) "water_param nueff Sedimentation:", nuice_sed
644
645	write(,) "Sed effective variance for CO2 clouds ?"
646	nuiceco2_sed=0.45
647	call getin_p("nuiceco2_sed",nuiceco2_sed)
648	write(,) "CO2 nueff Sedimentation:", nuiceco2_sed
649
650	write(,) "REF effective variance for CO2 clouds ?"
651	nuiceco2_ref=0.45
652	call getin_p("nuiceco2_ref",nuiceco2_ref)
653	write(,) "CO2 nueff Sedimentation:", nuiceco2_ref
654
655	write(,) "REF effective variance for water clouds ?"
656	nuice_ref=0.45
657	call getin_p("nuice_ref",nuice_ref)
658	write(,) "CO2 nueff Sedimentation:", nuice_ref
659
660
661	! ccn factor if no scavenging
662	write(,) "water param CCN reduc. factor ?"
663	ccn_factor = 4.5
664	call getin_p("ccn_factor",ccn_factor)
665	write(,)" ccn_factor = ",ccn_factor
666	write(,)"Careful: only used when microphys=F, otherwise"
667	write(,)"the contact parameter is used instead;"
668
669	! microphys
670	write(,)"Microphysical scheme for water-ice clouds?"
671	microphys=.false. ! default value
672	call getin_p("microphys",microphys)
673	write(,)" microphys = ",microphys
674
675	! supersat
676	write(,)"Allow super-saturation of water vapor?"
677	supersat=.true. ! default value
678	call getin_p("supersat",supersat)
679	write(,)"supersat = ",supersat
680
681	! microphysical parameter contact
682	write(,) "water contact parameter ?"
683	mteta = 0.95
684	call getin_p("mteta",mteta)
685	write(,) "mteta = ", mteta
686
687	! scavenging
688	write(,)"Dust scavenging by H2O/CO2 snowfall ?"
689	scavenging=.false. ! default value
690	call getin_p("scavenging",scavenging)
691	write(,)" scavenging = ",scavenging
692
693
694	! Test of incompatibility:
695	! if scavenging is used, then dustbin should be > 0
696
697	if ((microphys.and..not.doubleq).or.
698	& (microphys.and..not.water)) then
699	print*,'if microphys is used, then doubleq,'
700	print*,'and water must be used!'
701	call abort_physic(modname,
702	& "microphys requires water and doubleq",1)
703	endif
704	if (microphys.and..not.scavenging) then
705	print*,''
706	print*,'----------------WARNING-----------------'
707	print*,'microphys is used without scavenging !!!'
708	print*,'----------------WARNING-----------------'
709	print*,''
710	endif
711
712	if ((scavenging.and..not.microphys).or.
713	& (scavenging.and.(dustbin.lt.1)))then
714	print*,'if scavenging is used, then microphys'
715	print*,'must be used!'
716	call abort_physic(modname,
717	& "scavenging requires microphys",1)
718	endif
719
720	! Instantaneous scavenging by CO2
721	! -> expected to be replaced by scavenging with microphysics (flag scavenging) one day
722	write(,)"Dust scavenging by instantaneous CO2 snowfall ?"
723	scavco2cond=.false. ! default value
724	call getin_p("scavco2cond",scavco2cond)
725	write(,)" scavco2cond = ",scavco2cond
726	! Test of incompatibility:
727	! if scavco2cond is used, then dustbin should be > 0
728	if (scavco2cond.and.(dustbin.lt.1))then
729	print*,'if scavco2cond is used, then dustbin should be > 0'
730	call abort_physic(modname,
731	& "scavco2cond requires dustbin > 0",1)
732	endif
733	! if co2clouds is used, then there is no need for scavco2cond
734	if (co2clouds.and.scavco2cond) then
735	print*,''
736	print*,'----------------WARNING-----------------'
737	print*,' microphys scavenging is used so '
738	print*,' no need for scavco2cond !!! '
739	print*,'----------------WARNING-----------------'
740	print*,''
741	call abort_physic(modname,
742	& "incompatible co2cloud and scavco2cond options",1)
743	endif
744
745	! Test of incompatibility:
746
747	write(,) "Permanent water caps at poles ?",
748	& " .true. is RECOMMENDED"
749	write(,) "(with .true., North cap is a source of water ",
750	& "and South pole is a cold trap)"
751	caps=.true. ! default value
752	call getin_p("caps",caps)
753	write(,) " caps = ",caps
754
755	! albedo_h2o_ice
756	write(,) "water ice albedo ?"
757	albedo_h2o_ice=0.45
758	call getin_p("albedo_h2o_ice",albedo_h2o_ice)
759	write(,) " albedo_h2o_ice = ",albedo_h2o_ice
760	! inert_h2o_ice
761	write(,) "water ice thermal inertia ?"
762	inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
763	call getin_p("inert_h2o_ice",inert_h2o_ice)
764	write(,) " inert_h2o_ice = ",inert_h2o_ice
765	! frost_albedo_threshold
766	write(,) "frost thickness threshold for albedo ?"
767	frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2)
768	call getin_p("frost_albedo_threshold",
769	& frost_albedo_threshold)
770	write(,) " frost_albedo_threshold = ",
771	& frost_albedo_threshold
772
773	! call Titus crocus line -- DEFAULT IS NONE
774	write(,) "Titus crocus line ?"
775	tituscap=.false. ! default value
776	call getin_p("tituscap",tituscap)
777	write(,) "tituscap",tituscap
778
779	! Chemistry:
780	write(,) "photochemistry: include chemical species"
781	photochem=.false. ! default value
782	call getin_p("photochem",photochem)
783	write(,) " photochem = ",photochem
784
785	write(,) "Compute chemistry (if photochem is .true.)",
786	& "every ichemistry physics step (default: ichemistry=1)"
787	ichemistry=1
788	call getin_p("ichemistry",ichemistry)
789	write(,) " ichemistry = ",ichemistry
790
791
792	! SCATTERERS
793	write(,) "how many scatterers?"
794	naerkind=1 ! default value
795	call getin_p("naerkind",naerkind)
796	write(,)" naerkind = ",naerkind
797
798	! Test of incompatibility
799	c Logical tests for radiatively active water-ice clouds:
800	IF ( (activice.AND.(.NOT.water)).OR.
801	& (activice.AND.(naerkind.LT.2)) ) THEN
802	WRITE(,) 'If activice is TRUE, water has to be set'
803	WRITE(,) 'to TRUE, and "naerkind" must be at least'
804	WRITE(,) 'equal to 2.'
805	call abort_physic(modname,
806	& "radiatively active dust and water"//
807	& " require naerkind > 1",1)
808	ENDIF
809
810	!------------------------------------------
811	!------------------------------------------
812	! once naerkind is known allocate arrays
813	! -- we do it here and not in phys_var_init
814	! -- because we need to know naerkind first
815	CALL ini_scatterers(ngrid,nlayer)
816	!------------------------------------------
817	!------------------------------------------
818
819
820	c Please name the different scatterers here ----------------
821	name_iaer(1) = "dust_conrath" !! default choice is good old Conrath profile
822	IF (doubleq.AND.active) name_iaer(1) = "dust_doubleq" !! two-moment scheme
823	if (nq.gt.1) then
824	! trick to avoid problems compiling with 1 tracer
825	! and picky compilers who know name_iaer(2) is out of bounds
826	j=2
827	IF (rdstorm.AND..NOT.activice.AND..NOT.slpwind) name_iaer(2) =
828	& "stormdust_doubleq" !! storm dust two-moment scheme
829	IF (rdstorm.AND.water.AND.activice.AND..NOT.slpwind)
830	& name_iaer(3) = "stormdust_doubleq"
831	IF (slpwind.AND..NOT.activice.AND..NOT.rdstorm) name_iaer(2) =
832	& "topdust_doubleq" !! storm dust two-moment scheme
833	IF (slpwind.AND.water.AND.activice.AND..NOT.rdstorm)
834	& name_iaer(3) = "topdust_doubleq"
835	IF (rdstorm.AND.slpwind.AND..NOT.activice) THEN
836	name_iaer(2) = "stormdust_doubleq"
837	name_iaer(3) = "topdust_doubleq"
838	ENDIF
839	IF (rdstorm.AND.slpwind.AND.water.AND.activice) THEN
840	name_iaer(3) = "stormdust_doubleq"
841	name_iaer(4) = "topdust_doubleq"
842	ENDIF
843	IF (water.AND.activice) name_iaer(j) = "h2o_ice" !! radiatively-active clouds
844	IF (submicron.AND.active) name_iaer(j) = "dust_submicron" !! JBM experimental stuff
845	endif ! of if (nq.gt.1)
846
847	c ----------------------------------------------------------
848
849	! THERMOSPHERE
850
851	write(,) "call thermosphere ?"
852	callthermos=.false. ! default value
853	call getin_p("callthermos",callthermos)
854	write(,) " callthermos = ",callthermos
855
856
857	write(,) " water included without cycle ",
858	& "(only if water=.false.)"
859	thermoswater=.false. ! default value
860	call getin_p("thermoswater",thermoswater)
861	write(,) " thermoswater = ",thermoswater
862
863	write(,) "call thermal conduction ?",
864	& " (only if callthermos=.true.)"
865	callconduct=.false. ! default value
866	call getin_p("callconduct",callconduct)
867	write(,) " callconduct = ",callconduct
868
869	write(,) "call EUV heating ?",
870	& " (only if callthermos=.true.)"
871	calleuv=.false. ! default value
872	call getin_p("calleuv",calleuv)
873	write(,) " calleuv = ",calleuv
874
875	write(,) "call molecular viscosity ?",
876	& " (only if callthermos=.true.)"
877	callmolvis=.false. ! default value
878	call getin_p("callmolvis",callmolvis)
879	write(,) " callmolvis = ",callmolvis
880
881	write(,) "call molecular diffusion ?",
882	& " (only if callthermos=.true.)"
883	callmoldiff=.false. ! default value
884	call getin_p("callmoldiff",callmoldiff)
885	write(,) " callmoldiff = ",callmoldiff
886
887
888	write(,) "call thermospheric photochemistry ?",
889	& " (only if callthermos=.true.)"
890	thermochem=.false. ! default value
891	call getin_p("thermochem",thermochem)
892	write(,) " thermochem = ",thermochem
893
894	write(,) "Method to include solar variability"
895	write(,) "0-> fixed value of E10.7 (fixed_euv_value); ",
896	& "1-> daily evolution of E10.7 (for given solvaryear)"
897	solvarmod=1
898	call getin_p("solvarmod",solvarmod)
899	write(,) " solvarmod = ",solvarmod
900
901	write(,) "Fixed euv (for solvarmod==0) 10.7 value?"
902	write(,) " (min=80 , ave=140, max=320)"
903	fixed_euv_value=140 ! default value
904	call getin_p("fixed_euv_value",fixed_euv_value)
905	write(,) " fixed_euv_value = ",fixed_euv_value
906
907	write(,) "Solar variability as observed for MY: "
908	write(,) "Only if solvarmod=1"
909	solvaryear=24
910	call getin_p("solvaryear",solvaryear)
911	write(,) " solvaryear = ",solvaryear
912
913	write(,) "UV heating efficiency:",
914	& "measured values between 0.19 and 0.23 (Fox et al. 1996)",
915	& "lower values may be used to compensate low 15 um cooling"
916	euveff=0.21 !default value
917	call getin_p("euveff",euveff)
918	write(,) " euveff = ", euveff
919
920
921	if (.not.callthermos) then
922	if (thermoswater) then
923	print*,'if thermoswater is set, callthermos must be true'
924	call abort_physic(modname,
925	& "thermoswater requires callthermos",1)
926	endif
927	if (callconduct) then
928	print*,'if callconduct is set, callthermos must be true'
929	call abort_physic(modname,
930	& "callconduct requires callthermos",1)
931	endif
932	if (calleuv) then
933	print*,'if calleuv is set, callthermos must be true'
934	call abort_physic(modname,
935	& "calleuv requires callthermos",1)
936	endif
937	if (callmolvis) then
938	print*,'if callmolvis is set, callthermos must be true'
939	call abort_physic(modname,
940	& "callmolvis requires callthermos",1)
941	endif
942	if (callmoldiff) then
943	print*,'if callmoldiff is set, callthermos must be true'
944	call abort_physic(modname,
945	& "callmoldiff requires callthermos",1)
946	endif
947	if (thermochem) then
948	print*,'if thermochem is set, callthermos must be true'
949	call abort_physic(modname,
950	& "thermochem requires callthermos",1)
951	endif
952	endif
953
954	! Test of incompatibility:
955	! if photochem is used, then water should be used too
956
957	if (photochem.and..not.water) then
958	print*,'if photochem is used, water should be used too'
959	call abort_physic(modname,
960	& "photochem requires water",1)
961	endif
962
963	! if callthermos is used, then thermoswater should be used too
964	! (if water not used already)
965
966	if (callthermos .and. .not.water) then
967	if (callthermos .and. .not.thermoswater) then
968	print*,'if callthermos is used, water or thermoswater
969	& should be used too'
970	call abort_physic(modname,
971	& "callthermos requires water or thermoswater",1)
972	endif
973	endif
974
975	PRINT*,'--------------------------------------------'
976	PRINT*
977	PRINT*
978	ELSE
979	write(,)
980	write(,) 'Cannot read file callphys.def. Is it here ?'
981	call abort_physic(modname,
982	& "missing callphys.def file",1)
983	ENDIF
984
985	8000 FORMAT(t5,a12,l8)
986	8001 FORMAT(t5,a12,i8)
987
988	PRINT*
989	PRINT*,'conf_phys: daysec',daysec
990	PRINT*
991	PRINT*,'conf_phys: The radiative transfer is computed:'
992	PRINT*,' each ',iradia,' physical time-step'
993	PRINT,' or each ',iradiadtphys,' seconds'
994	PRINT*
995	! --------------------------------------------------------------
996	! Managing the Longwave radiative transfer
997	! --------------------------------------------------------------
998
999	! In most cases, the run just use the following values :
1000	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1001	callemis=.true.
1002	! ilwd=10*int(daysec/dtphys) ! bug before 22/10/01
1003	ilwd=1
1004	ilwn=1 !2
1005	ilwb=1 !2
1006	linear=.true.
1007	ncouche=3
1008	alphan=0.4
1009	semi=0
1010
1011	! BUT people working hard on the LW may want to read them in 'radia.def'
1012	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1013
1014	OPEN(99,file='radia.def',status='old',form='formatted'
1015	. ,iostat=ierr)
1016	IF(ierr.EQ.0) THEN
1017	write(,) 'conf_phys: Reading radia.def !!!'
1018	READ(99,fmt='(a)') ch1
1019	READ(99,*) callemis
1020	WRITE(*,8000) ch1,callemis
1021
1022	READ(99,fmt='(a)') ch1
1023	READ(99,*) iradia
1024	WRITE(*,8001) ch1,iradia
1025
1026	READ(99,fmt='(a)') ch1
1027	READ(99,*) ilwd
1028	WRITE(*,8001) ch1,ilwd
1029
1030	READ(99,fmt='(a)') ch1
1031	READ(99,*) ilwn
1032	WRITE(*,8001) ch1,ilwn
1033
1034	READ(99,fmt='(a)') ch1
1035	READ(99,*) linear
1036	WRITE(*,8000) ch1,linear
1037
1038	READ(99,fmt='(a)') ch1
1039	READ(99,*) ncouche
1040	WRITE(*,8001) ch1,ncouche
1041
1042	READ(99,fmt='(a)') ch1
1043	READ(99,*) alphan
1044	WRITE(,) ch1,alphan
1045
1046	READ(99,fmt='(a)') ch1
1047	READ(99,*) ilwb
1048	WRITE(*,8001) ch1,ilwb
1049
1050
1051	READ(99,fmt='(a)') ch1
1052	READ(99,'(l1)') callg2d
1053	WRITE(*,8000) ch1,callg2d
1054
1055	READ(99,fmt='(a)') ch1
1056	READ(99,*) semi
1057	WRITE(,) ch1,semi
1058	end if
1059	CLOSE(99)
1060
1061	END

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

Download in other formats:

Original Format