Context Navigation

conf_phys.F @ 2156

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

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

Download in other formats:

Original Format