Context Navigation

conf_phys.F @ 2947

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

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

Download in other formats:

Original Format