Context Navigation

conf_phys.F @ 3111

Last change on this file since 3111 was 3111, checked in by llange, 13 months ago
MARS PCM 1) Following r-3098, adding d_coeff in the OMPthreadprivate 2) Introduce a temporary "old_wsublimation_scheme": when set to true (by default), the water frost sublimation is computed as usual. Else, it is computed as rho Ch U (q-qsat) and not rho Cd U (q-qsat) LL
File size: 46.0 KB

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

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

Download in other formats:

Original Format