Context Navigation

conf_phys.F @ 2740

Last change on this file since 2740 was 2661, checked in by emillour, 3 years ago
Mars GCM: Follow-up to removal of nuice_ref default value in initracer. Put nuice_ref=01 as default value in conf_phys to maintain default behaviour. EM
File size: 44.7 KB

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

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

Download in other formats:

Original Format