Context Navigation

conf_phys.F @ 2740

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

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

Download in other formats:

Original Format