Context Navigation

conf_phys.F @ 2199

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

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

Download in other formats:

Original Format