Context Navigation

conf_phys.F @ 3027

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

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

Download in other formats:

Original Format