Context Navigation

conf_phys.F @ 2627

Last change on this file since 2627 was 2627, checked in by romain.vande, 3 years ago
LMDZ_MARS: Small OpenMP fixes in conf_phys for reading radia.def with ifort RV
File size: 43.5 KB

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

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

Download in other formats:

Original Format