Context Navigation

conf_phys.F @ 3712

Last change on this file since 3712 was 3712, checked in by emillour, 3 months ago
Mars PCM: Add extra checking in nirco2abs to ensure having nircorr==1 when o and co2 tracers are available. EM
File size: 50.3 KB

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

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

Download in other formats:

Original Format