Context Navigation

conf_phys.F @ 3026

Last change on this file since 3026 was 3008, checked in by emillour, 18 months ago
Mars PCM: Some code cleanup around microphysics. Turn microphys.h into module microphys_h.F90, and while at it also turn nuclea.F, growthrate.F90 and massflowrateco2.F90 into modules. EM
File size: 45.4 KB

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

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

Download in other formats:

Original Format