Context Navigation

conf_phys.F @ 2932

Last change on this file since 2932 was 2916, checked in by emillour, 22 months ago
Mars PCM Add a "diagsoil" flag to trigger outputing a diagsoil.nc file (default is diagsoil=.false.) EM
File size: 44.5 KB

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

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

Download in other formats:

Original Format