Context Navigation

conf_phys.F @ 2179

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

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

Download in other formats:

Original Format