Context Navigation

conf_phys.F @ 2538

Last change on this file since 2538 was 2522, checked in by jnaar, 4 years ago
Added flag for temperature-dependant water contact parameter. Default is false for now ! JN
File size: 42.2 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 ! default value
709	temp_dependant_m = .false. ! default value
710	call getin_p("temp_dependant_m",temp_dependant_m)
711	if (temp_dependant_m) then
712	print*,'You have chosen a temperature-dependant water'
713	print*,'contact parameter ! From Maattanen et al. 2014'
714	else if (.not.temp_dependant_m) then
715	print*,'Water contact parameter is constant'
716	call getin_p("mteta",mteta)
717	write(,) "mteta = ", mteta
718	endif
719
720	! scavenging
721	write(,)"Dust scavenging by H2O/CO2 snowfall ?"
722	scavenging=.false. ! default value
723	call getin_p("scavenging",scavenging)
724	write(,)" scavenging = ",scavenging
725
726
727	! Test of incompatibility:
728	! if scavenging is used, then dustbin should be > 0
729
730	if ((microphys.and..not.doubleq).or.
731	& (microphys.and..not.water)) then
732	print*,'if microphys is used, then doubleq,'
733	print*,'and water must be used!'
734	call abort_physic(modname,
735	& "microphys requires water and doubleq",1)
736	endif
737	if (microphys.and..not.scavenging) then
738	print*,''
739	print*,'----------------WARNING-----------------'
740	print*,'microphys is used without scavenging !!!'
741	print*,'----------------WARNING-----------------'
742	print*,''
743	endif
744
745	if ((scavenging.and..not.microphys).or.
746	& (scavenging.and.(dustbin.lt.1)))then
747	print*,'if scavenging is used, then microphys'
748	print*,'must be used!'
749	call abort_physic(modname,
750	& "scavenging requires microphys",1)
751	endif
752
753	! Instantaneous scavenging by CO2
754	! -> expected to be replaced by scavenging with microphysics (flag scavenging) one day
755	write(,)"Dust scavenging by instantaneous CO2 snowfall ?"
756	scavco2cond=.false. ! default value
757	call getin_p("scavco2cond",scavco2cond)
758	write(,)" scavco2cond = ",scavco2cond
759	! Test of incompatibility:
760	! if scavco2cond is used, then dustbin should be > 0
761	if (scavco2cond.and.(dustbin.lt.1))then
762	print*,'if scavco2cond is used, then dustbin should be > 0'
763	call abort_physic(modname,
764	& "scavco2cond requires dustbin > 0",1)
765	endif
766	! if co2clouds is used, then there is no need for scavco2cond
767	if (co2clouds.and.scavco2cond) then
768	print*,''
769	print*,'----------------WARNING-----------------'
770	print*,' microphys scavenging is used so '
771	print*,' no need for scavco2cond !!! '
772	print*,'----------------WARNING-----------------'
773	print*,''
774	call abort_physic(modname,
775	& "incompatible co2cloud and scavco2cond options",1)
776	endif
777
778	! Test of incompatibility:
779
780	write(,) "Permanent water caps at poles ?",
781	& " .true. is RECOMMENDED"
782	write(,) "(with .true., North cap is a source of water ",
783	& "and South pole is a cold trap)"
784	caps=.true. ! default value
785	call getin_p("caps",caps)
786	write(,) " caps = ",caps
787
788	! JN : now separated between albedo_h2o_cap and
789	! albedo_h2o_frost. Retrocompatible with old
790	! callphys.def with albedo_h2o_ice
791	write(,) "water ice albedo ? Old settings use ",
792	& "albedo_h2o_ice, new settings use ",
793	& "albedo_h2o_cap and albedo_h2o_frost "
794	albedo_h2o_cap=0.35
795	albedo_h2o_frost=0.35
796	call getin_p("albedo_h2o_ice",albedo_h2o_cap)
797	albedo_h2o_frost=albedo_h2o_cap
798	call getin_p("albedo_h2o_cap",albedo_h2o_cap)
799	write(,) " albedo_h2o_cap = ",albedo_h2o_cap
800	call getin_p("albedo_h2o_frost",albedo_h2o_frost)
801	write(,) " albedo_h2o_frost = ",albedo_h2o_frost
802
803	! Northern polar cap albedo (JN 2021)
804	write(,)"Watercaptag albedo is unchanged by water frost",
805	& " deposition (default is false)"
806	cap_albedo=.false. ! default value
807	call getin_p("cap_albedo",cap_albedo)
808	write(,)"cap_albedo = ",cap_albedo
809
810	! inert_h2o_ice
811	write(,) "water ice thermal inertia ?"
812	inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
813	call getin_p("inert_h2o_ice",inert_h2o_ice)
814	write(,) " inert_h2o_ice = ",inert_h2o_ice
815	! frost_albedo_threshold
816	write(,) "frost thickness threshold for albedo ?"
817	frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2)
818	call getin_p("frost_albedo_threshold",
819	& frost_albedo_threshold)
820	write(,) " frost_albedo_threshold = ",
821	& frost_albedo_threshold
822
823	! call Titus crocus line -- DEFAULT IS NONE
824	write(,) "Titus crocus line ?"
825	tituscap=.false. ! default value
826	call getin_p("tituscap",tituscap)
827	write(,) "tituscap",tituscap
828
829	! Chemistry:
830	write(,) "photochemistry: include chemical species"
831	photochem=.false. ! default value
832	call getin_p("photochem",photochem)
833	write(,) " photochem = ",photochem
834
835	write(,) "Compute chemistry (if photochem is .true.)",
836	& "every ichemistry physics step (default: ichemistry=1)"
837	ichemistry=1
838	call getin_p("ichemistry",ichemistry)
839	write(,) " ichemistry = ",ichemistry
840
841
842	! SCATTERERS
843	write(,) "how many scatterers?"
844	naerkind=1 ! default value
845	call getin_p("naerkind",naerkind)
846	write(,)" naerkind = ",naerkind
847
848	! Test of incompatibility
849	c Logical tests for radiatively active water-ice clouds:
850	IF ( (activice.AND.(.NOT.water)).OR.
851	& (activice.AND.(naerkind.LT.2)) ) THEN
852	WRITE(,) 'If activice is TRUE, water has to be set'
853	WRITE(,) 'to TRUE, and "naerkind" must be at least'
854	WRITE(,) 'equal to 2.'
855	call abort_physic(modname,
856	& "radiatively active dust and water"//
857	& " require naerkind > 1",1)
858	ENDIF
859
860	!------------------------------------------
861	!------------------------------------------
862	! once naerkind is known allocate arrays
863	! -- we do it here and not in phys_var_init
864	! -- because we need to know naerkind first
865	CALL ini_scatterers(ngrid,nlayer)
866	!------------------------------------------
867	!------------------------------------------
868
869
870	c Please name the different scatterers here ----------------
871	name_iaer(1) = "dust_conrath" !! default choice is good old Conrath profile
872	IF (doubleq.AND.active) name_iaer(1) = "dust_doubleq" !! two-moment scheme
873
874	if (nq.gt.1) then
875	! trick to avoid problems compiling with 1 tracer
876	! and picky compilers who know name_iaer(2) is out of bounds
877	j=2
878	IF (rdstorm.AND..NOT.activice.AND..NOT.slpwind) then
879	name_iaer(j) = "stormdust_doubleq" !! storm dust two-moment scheme
880	j = j+1
881	END IF
882
883	IF (rdstorm.AND.water.AND.activice.AND..NOT.slpwind) then
884	name_iaer(j) = "stormdust_doubleq"
885	j = j+1
886	END IF
887
888	IF (slpwind.AND..NOT.activice.AND..NOT.rdstorm) then
889	name_iaer(j) = "topdust_doubleq" !! storm dust two-moment scheme
890	j = j+1
891	END IF
892
893	IF (slpwind.AND.water.AND.activice.AND..NOT.rdstorm) then
894	name_iaer(j) = "topdust_doubleq"
895	j = j+1
896	END IF
897
898	IF (rdstorm.AND.slpwind.AND..NOT.activice) THEN
899	name_iaer(j) = "stormdust_doubleq"
900	name_iaer(j+1) = "topdust_doubleq"
901	j = j+2
902	ENDIF
903
904	IF (rdstorm.AND.slpwind.AND.water.AND.activice) THEN
905	name_iaer(j) = "stormdust_doubleq"
906	name_iaer(j+1) = "topdust_doubleq"
907	j = j+2
908	ENDIF
909
910	IF (water.AND.activice) then
911	name_iaer(j) = "h2o_ice" !! radiatively-active clouds
912	j = j+1
913	END IF
914
915	IF (co2clouds.AND.activeco2ice) then
916	name_iaer(j) = "co2_ice" !! radiatively-active co2 clouds
917	j = j+1
918	ENDIF
919
920	IF (submicron.AND.active) then
921	name_iaer(j) = "dust_submicron" !! JBM experimental stuff
922	j = j+1
923	ENDIF
924	endif ! of if (nq.gt.1)
925	c ----------------------------------------------------------
926
927	! THERMOSPHERE
928
929	write(,) "call thermosphere ?"
930	callthermos=.false. ! default value
931	call getin_p("callthermos",callthermos)
932	write(,) " callthermos = ",callthermos
933
934
935	write(,) " water included without cycle ",
936	& "(only if water=.false.)"
937	thermoswater=.false. ! default value
938	call getin_p("thermoswater",thermoswater)
939	write(,) " thermoswater = ",thermoswater
940
941	write(,) "call thermal conduction ?",
942	& " (only if callthermos=.true.)"
943	callconduct=.false. ! default value
944	call getin_p("callconduct",callconduct)
945	write(,) " callconduct = ",callconduct
946
947	write(,) "call EUV heating ?",
948	& " (only if callthermos=.true.)"
949	calleuv=.false. ! default value
950	call getin_p("calleuv",calleuv)
951	write(,) " calleuv = ",calleuv
952
953	write(,) "call molecular viscosity ?",
954	& " (only if callthermos=.true.)"
955	callmolvis=.false. ! default value
956	call getin_p("callmolvis",callmolvis)
957	write(,) " callmolvis = ",callmolvis
958
959	write(,) "call molecular diffusion ?",
960	& " (only if callthermos=.true.)"
961	callmoldiff=.false. ! default value
962	call getin_p("callmoldiff",callmoldiff)
963	write(,) " callmoldiff = ",callmoldiff
964
965
966	write(,) "call thermospheric photochemistry ?",
967	& " (only if callthermos=.true.)"
968	thermochem=.false. ! default value
969	call getin_p("thermochem",thermochem)
970	write(,) " thermochem = ",thermochem
971
972	write(,) "Method to include solar variability"
973	write(,) "0-> fixed value of E10.7 (fixed_euv_value); ",
974	& "1-> daily evolution of E10.7 (for given solvaryear)"
975	solvarmod=1
976	call getin_p("solvarmod",solvarmod)
977	write(,) " solvarmod = ",solvarmod
978
979	write(,) "Fixed euv (for solvarmod==0) 10.7 value?"
980	write(,) " (min=80 , ave=140, max=320)"
981	fixed_euv_value=140 ! default value
982	call getin_p("fixed_euv_value",fixed_euv_value)
983	write(,) " fixed_euv_value = ",fixed_euv_value
984
985	write(,) "Solar variability as observed for MY: "
986	write(,) "Only if solvarmod=1"
987	solvaryear=24
988	call getin_p("solvaryear",solvaryear)
989	write(,) " solvaryear = ",solvaryear
990
991	write(,) "UV heating efficiency:",
992	& "measured values between 0.19 and 0.23 (Fox et al. 1996)",
993	& "lower values may be used to compensate low 15 um cooling"
994	euveff=0.21 !default value
995	call getin_p("euveff",euveff)
996	write(,) " euveff = ", euveff
997
998
999	if (.not.callthermos) then
1000	if (thermoswater) then
1001	print*,'if thermoswater is set, callthermos must be true'
1002	call abort_physic(modname,
1003	& "thermoswater requires callthermos",1)
1004	endif
1005	if (callconduct) then
1006	print*,'if callconduct is set, callthermos must be true'
1007	call abort_physic(modname,
1008	& "callconduct requires callthermos",1)
1009	endif
1010	if (calleuv) then
1011	print*,'if calleuv is set, callthermos must be true'
1012	call abort_physic(modname,
1013	& "calleuv requires callthermos",1)
1014	endif
1015	if (callmolvis) then
1016	print*,'if callmolvis is set, callthermos must be true'
1017	call abort_physic(modname,
1018	& "callmolvis requires callthermos",1)
1019	endif
1020	if (callmoldiff) then
1021	print*,'if callmoldiff is set, callthermos must be true'
1022	call abort_physic(modname,
1023	& "callmoldiff requires callthermos",1)
1024	endif
1025	if (thermochem) then
1026	print*,'if thermochem is set, callthermos must be true'
1027	call abort_physic(modname,
1028	& "thermochem requires callthermos",1)
1029	endif
1030	endif
1031
1032	! Test of incompatibility:
1033	! if photochem is used, then water should be used too
1034
1035	if (photochem.and..not.water) then
1036	print*,'if photochem is used, water should be used too'
1037	call abort_physic(modname,
1038	& "photochem requires water",1)
1039	endif
1040
1041	! if callthermos is used, then thermoswater should be used too
1042	! (if water not used already)
1043
1044	if (callthermos .and. .not.water) then
1045	if (callthermos .and. .not.thermoswater) then
1046	print*,'if callthermos is used, water or thermoswater
1047	& should be used too'
1048	call abort_physic(modname,
1049	& "callthermos requires water or thermoswater",1)
1050	endif
1051	endif
1052
1053	PRINT*,'--------------------------------------------'
1054	PRINT*
1055	PRINT*
1056	ELSE
1057	write(,)
1058	write(,) 'Cannot read file callphys.def. Is it here ?'
1059	call abort_physic(modname,
1060	& "missing callphys.def file",1)
1061	ENDIF
1062
1063	8000 FORMAT(t5,a12,l8)
1064	8001 FORMAT(t5,a12,i8)
1065
1066	PRINT*
1067	PRINT*,'conf_phys: daysec',daysec
1068	PRINT*
1069	PRINT*,'conf_phys: The radiative transfer is computed:'
1070	PRINT*,' each ',iradia,' physical time-step'
1071	PRINT,' or each ',iradiadtphys,' seconds'
1072	PRINT*
1073	! --------------------------------------------------------------
1074	! Managing the Longwave radiative transfer
1075	! --------------------------------------------------------------
1076
1077	! In most cases, the run just use the following values :
1078	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1079	callemis=.true.
1080	! ilwd=10*int(daysec/dtphys) ! bug before 22/10/01
1081	ilwd=1
1082	ilwn=1 !2
1083	ilwb=1 !2
1084	linear=.true.
1085	ncouche=3
1086	alphan=0.4
1087	semi=0
1088
1089	! BUT people working hard on the LW may want to read them in 'radia.def'
1090	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1091
1092	OPEN(99,file='radia.def',status='old',form='formatted'
1093	. ,iostat=ierr)
1094	IF(ierr.EQ.0) THEN
1095	write(,) 'conf_phys: Reading radia.def !!!'
1096	READ(99,fmt='(a)') ch1
1097	READ(99,*) callemis
1098	WRITE(*,8000) ch1,callemis
1099
1100	READ(99,fmt='(a)') ch1
1101	READ(99,*) iradia
1102	WRITE(*,8001) ch1,iradia
1103
1104	READ(99,fmt='(a)') ch1
1105	READ(99,*) ilwd
1106	WRITE(*,8001) ch1,ilwd
1107
1108	READ(99,fmt='(a)') ch1
1109	READ(99,*) ilwn
1110	WRITE(*,8001) ch1,ilwn
1111
1112	READ(99,fmt='(a)') ch1
1113	READ(99,*) linear
1114	WRITE(*,8000) ch1,linear
1115
1116	READ(99,fmt='(a)') ch1
1117	READ(99,*) ncouche
1118	WRITE(*,8001) ch1,ncouche
1119
1120	READ(99,fmt='(a)') ch1
1121	READ(99,*) alphan
1122	WRITE(,) ch1,alphan
1123
1124	READ(99,fmt='(a)') ch1
1125	READ(99,*) ilwb
1126	WRITE(*,8001) ch1,ilwb
1127
1128
1129	READ(99,fmt='(a)') ch1
1130	READ(99,'(l1)') callg2d
1131	WRITE(*,8000) ch1,callg2d
1132
1133	READ(99,fmt='(a)') ch1
1134	READ(99,*) semi
1135	WRITE(,) ch1,semi
1136	end if
1137	CLOSE(99)
1138
1139	END

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

Download in other formats:

Original Format