Context Navigation

conf_phys.F @ 2341

Last change on this file since 2341 was 2312, checked in by lrossi, 5 years ago

MARS GCM:
Implementation of HDO in the physics
New tracers hdo_vap and hdo_ice are added. Cf. traceur.def.isotopes in deftank for exemple of traceur.def.
All HDO related computations are activated by flag hdo=.true. in callphys.def. (see callphys.def.hdo in deftank for an example)
Additional option hdofrac (true by default) allows for turning on/off fractionation (for tests).
For now, only functional with simplified cloud scheme (so microphys=.false. and activice=.false. also recommended).
Initialisation of start files can be done with option 'inihdo' in newstart.
LR

File size: 39.5 KB

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

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

Download in other formats:

Original Format