Context Navigation

conf_phys.F @ 2879

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

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

Download in other formats:

Original Format