Context Navigation

conf_phys.F @ 2559

Last change on this file since 2559 was 2559, checked in by emillour, 3 years ago

Mars GCM:
Some code cleanup and refactoring around wstats:

turn wstats.F90 in a module
remove no useless statto_mod.F90
incorporate auxilliary inistats and mkstats routines in wstats_mod.F90
move flag "callstats" from callkeys.h to being a module variable of wstats_mod

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

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

Download in other formats:

Original Format