Context Navigation

conf_phys.F @ 2997

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

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

Download in other formats:

Original Format