Context Navigation

conf_phys.F @ 2825

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

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

Download in other formats:

Original Format