Context Navigation

conf_phys.F @ 1962

Last change on this file since 1962 was 1918, checked in by emillour, 7 years ago

Mars GCM:
Code cleanup:

remove "comorbit.h" since it is no longer used.
turn "datafile.h" into module datafile_mod.F90 (and rename variable "datafile" as "datadir" since it stores the path to the datafile directory).

File size: 30.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
	282
	283	write(,)"rad.transfer is computed every iradia",
	284	& " physical timestep"
	285	iradia=1 ! default value
	286	call getin("iradia",iradia)
	287	write(,)" iradia = ",iradia
	288
	289
	290	write(,)"Output of the exchange coefficient mattrix ?",
	291	& "(for diagnostics only)"
	292	callg2d=.false. ! default value
	293	call getin("callg2d",callg2d)
	294	write(,)" callg2d = ",callg2d
	295
	296	write(,)"Rayleigh scattering : (should be .false. for now)"
	297	rayleigh=.false.
	298	call getin("rayleigh",rayleigh)
	299	write(,)" rayleigh = ",rayleigh
	300
	301
	302	! TRACERS:
	303
[358]	304	! dustbin
[42]	305	write(,)"Transported dust ? (if >0, use 'dustbin' dust bins)"
	306	dustbin=0 ! default value
	307	call getin("dustbin",dustbin)
	308	write(,)" dustbin = ",dustbin
[358]	309	! active
[42]	310	write(,)"Radiatively active dust ? (matters if dustbin>0)"
	311	active=.false. ! default value
	312	call getin("active",active)
	313	write(,)" active = ",active
	314
	315	! Test of incompatibility:
	316	! if active is used, then dustbin should be > 0
	317
	318	if (active.and.(dustbin.lt.1)) then
	319	print*,'if active is used, then dustbin should > 0'
	320	stop
	321	endif
[358]	322	! doubleq
[42]	323	write(,)"use mass and number mixing ratios to predict",
	324	& " dust size ?"
	325	doubleq=.false. ! default value
	326	call getin("doubleq",doubleq)
	327	write(,)" doubleq = ",doubleq
[358]	328	! submicron
[42]	329	submicron=.false. ! default value
	330	call getin("submicron",submicron)
	331	write(,)" submicron = ",submicron
	332
	333	! Test of incompatibility:
	334	! if doubleq is used, then dustbin should be 2
	335
	336	if (doubleq.and.(dustbin.ne.2)) then
	337	print*,'if doubleq is used, then dustbin should be 2'
	338	stop
	339	endif
[1036]	340	if (doubleq.and.submicron.and.(nq.LT.3)) then
[42]	341	print*,'If doubleq is used with a submicron tracer,'
	342	print*,' then the number of tracers has to be'
	343	print*,' larger than 3.'
	344	stop
	345	endif
[1088]	346
[358]	347	! lifting
[42]	348	write(,)"dust lifted by GCM surface winds ?"
	349	lifting=.false. ! default value
	350	call getin("lifting",lifting)
	351	write(,)" lifting = ",lifting
	352
	353	! Test of incompatibility:
	354	! if lifting is used, then dustbin should be > 0
	355
	356	if (lifting.and.(dustbin.lt.1)) then
	357	print*,'if lifting is used, then dustbin should > 0'
	358	stop
	359	endif
[1088]	360
	361	! free evolving dust
	362	! freedust=true just says that there is no lifting and no dust opacity scaling.
	363	write(,)"dust lifted by GCM surface winds ?"
	364	freedust=.false. ! default value
	365	call getin("freedust",freedust)
	366	write(,)" freedust = ",freedust
	367	if (freedust.and..not.doubleq) then
	368	print*,'freedust should be used with doubleq !'
	369	stop
	370	endif
[1264]	371	#ifndef MESOSCALE
	372	! this test is valid in GCM case
	373	! ... not in mesoscale case, for we want to activate mesoscale lifting
[1088]	374	if (freedust.and.lifting) then
	375	print*,'if freedust is used, then lifting should not be used'
	376	print*,'lifting forced to false !!'
	377	lifting=.false.
	378	endif
[1264]	379	#endif
[1088]	380
[1353]	381	! Dust IR opacity
	382	write(,)" Wavelength for infrared opacity of dust ?"
	383	write(,)" Choices are:"
	384	write(,)" tes --- > 9.3 microns [default]"
	385	write(,)" mcs --- > 21.6 microns"
	386	!
	387	! WARNING WARNING WARNING WARNING WARNING WARNING
	388	!
	389	! BEFORE ADDING A NEW VALUE, BE SURE THAT THE
	390	! CORRESPONDING WAVELENGTH IS IN THE LOOKUP TABLE,
	391	! OR AT LEAST NO TO FAR, TO AVOID FALLACIOUS INTERPOLATIONS.
	392	!
	393	dustiropacity="tes" !default value - is expected to shift to mcs one day
	394	call getin("dustiropacity",dustiropacity)
	395	write(,)" dustiropacity = ",trim(dustiropacity)
	396	select case (trim(dustiropacity))
	397	case ("tes")
	398	dustrefir = 9.3E-6
	399	case ("mcs")
	400	dustrefir = 21.6E-6
	401	case default
	402	write(,) trim(dustiropacity),
	403	& " is not a valid option for dustiropacity"
	404	stop
	405	end select
	406
[358]	407	! callddevil
[42]	408	write(,)" dust lifted by dust devils ?"
	409	callddevil=.false. !default value
	410	call getin("callddevil",callddevil)
	411	write(,)" callddevil = ",callddevil
	412
	413	! Test of incompatibility:
	414	! if dustdevil is used, then dustbin should be > 0
	415
	416	if (callddevil.and.(dustbin.lt.1)) then
	417	print*,'if dustdevil is used, then dustbin should > 0'
	418	stop
	419	endif
[358]	420	! sedimentation
[42]	421	write(,) "Gravitationnal sedimentation ?"
	422	sedimentation=.true. ! default value
	423	call getin("sedimentation",sedimentation)
	424	write(,) " sedimentation = ",sedimentation
[358]	425	! activice
[42]	426	write(,) "Radiatively active transported atmospheric ",
	427	& "water ice ?"
	428	activice=.false. ! default value
	429	call getin("activice",activice)
	430	write(,) " activice = ",activice
[358]	431	! water
[42]	432	write(,) "Compute water cycle ?"
	433	water=.false. ! default value
	434	call getin("water",water)
	435	write(,) " water = ",water
	436
[1711]	437	! sub-grid cloud fraction: fixed
	438	write(,) "Fixed cloud fraction?"
	439	CLFfixval=1.0 ! default value
	440	call getin("CLFfixval",CLFfixval)
	441	write(,) "CLFfixval=",CLFfixval
	442	! sub-grid cloud fraction: varying
	443	write(,) "Use partial nebulosity?"
	444	CLFvarying=.false. ! default value
	445	call getin("CLFvarying",CLFvarying)
	446	write(,)"CLFvarying=",CLFvarying
	447
[1617]	448	!CO2 clouds scheme?
[1720]	449	write(,) "Compute CO2 clouds (implies microphysical scheme)?"
[1617]	450	co2clouds=.false. ! default value
	451	call getin("co2clouds",co2clouds)
	452	write(,) " co2clouds = ",co2clouds
[1720]	453	!Can water ice particles serve as CCN for CO2clouds
	454	write(,) "Use water ice as CO2 clouds CCN ?"
	455	co2useh2o=.false. ! default value
	456	call getin("co2useh2o",co2useh2o)
	457	write(,) " co2useh2o = ",co2useh2o
	458	!Do we allow a supply of meteoritic paricles to serve as CO2 ice CCN?
	459	write(,) "Supply meteoritic particle for CO2 clouds ?"
	460	meteo_flux=.false. !Default value
	461	call getin("meteo_flux",meteo_flux)
	462	write(,) " meteo_flux = ",meteo_flux
	463	!Do we allow a sub-grid temperature distribution for the CO2 microphysics
	464	write(,) "sub-grid temperature distribution for CO2 clouds?"
	465	CLFvaryingCO2=.false. !Default value
	466	call getin("CLFvaryingCO2",CLFvaryingCO2)
	467	write(,) " CLFvaryingCO2 = ",CLFvaryingCO2
	468	!Amplitude of the sub-grid temperature distribution for the CO2 microphysics
	469	write(,) "sub-grid temperature amplitude for CO2 clouds?"
	470	spantCO2=0 !Default value
	471	call getin("spantCO2",spantCO2)
	472	write(,) " spantCO2 = ",spantCO2
[1818]	473	!Do you want to filter the sub-grid T distribution by a Saturation index?
	474	write(,) "filter sub-grid temperature by Saturation index?"
	475	satindexco2=.true.
	476	call getin("satindexco2",satindexco2)
	477	write(,) " satindexco2 = ",satindexco2
[1720]	478
[1818]	479
[833]	480	! thermal inertia feedback
	481	write(,) "Activate the thermal inertia feedback ?"
	482	tifeedback=.false. ! default value
	483	call getin("tifeedback",tifeedback)
	484	write(,) " tifeedback = ",tifeedback
	485
[42]	486	! Test of incompatibility:
	487
[833]	488	if (tifeedback.and..not.water) then
	489	print*,'if tifeedback is used,'
	490	print*,'water should be used too'
	491	stop
	492	endif
	493
	494	if (tifeedback.and..not.callsoil) then
	495	print*,'if tifeedback is used,'
	496	print*,'callsoil should be used too'
	497	stop
	498	endif
	499
[42]	500	if (activice.and..not.water) then
	501	print*,'if activice is used, water should be used too'
	502	stop
	503	endif
	504
	505	if (water.and..not.tracer) then
	506	print*,'if water is used, tracer should be used too'
	507	stop
	508	endif
[420]	509
[455]	510	! water ice clouds effective variance distribution for sedimentaion
[1617]	511	write(,) "Sed effective variance for water ice clouds ?"
[455]	512	nuice_sed=0.45
	513	call getin("nuice_sed",nuice_sed)
	514	write(,) "water_param nueff Sedimentation:", nuice_sed
[1617]	515
	516	write(,) "Sed effective variance for CO2 clouds ?"
	517	nuiceco2_sed=0.45
	518	call getin("nuiceco2_sed",nuiceco2_sed)
	519	write(,) "CO2 nueff Sedimentation:", nuiceco2_sed
	520
	521	write(,) "REF effective variance for CO2 clouds ?"
	522	nuiceco2_ref=0.45
	523	call getin("nuiceco2_ref",nuiceco2_ref)
	524	write(,) "CO2 nueff Sedimentation:", nuiceco2_ref
	525
	526	write(,) "REF effective variance for water clouds ?"
	527	nuice_ref=0.45
	528	call getin("nuice_ref",nuice_ref)
	529	write(,) "CO2 nueff Sedimentation:", nuice_ref
	530
	531
[420]	532	! ccn factor if no scavenging
[1617]	533	write(,) "water param CCN reduc. factor ?"
[420]	534	ccn_factor = 4.5
	535	call getin("ccn_factor",ccn_factor)
	536	write(,)" ccn_factor = ",ccn_factor
	537	write(,)"Careful: only used when microphys=F, otherwise"
	538	write(,)"the contact parameter is used instead;"
	539
[1720]	540	! microphys
	541	write(,)"Microphysical scheme for water-ice clouds?"
	542	microphys=.false. ! default value
	543	call getin("microphys",microphys)
	544	write(,)" microphys = ",microphys
[1655]	545
[1720]	546	! supersat
	547	write(,)"Allow super-saturation of water vapor?"
	548	supersat=.true. ! default value
	549	call getin("supersat",supersat)
	550	write(,)"supersat = ",supersat
[42]	551
[358]	552	! microphysical parameter contact
[1720]	553	write(,) "water contact parameter ?"
	554	mteta = 0.95
	555	call getin("mteta",mteta)
	556	write(,) "mteta = ", mteta
	557
[358]	558	! scavenging
[1720]	559	write(,)"Dust scavenging by H2O/CO2 snowfall ?"
	560	scavenging=.false. ! default value
	561	call getin("scavenging",scavenging)
	562	write(,)" scavenging = ",scavenging
[358]	563
	564
[42]	565	! Test of incompatibility:
[358]	566	! if scavenging is used, then dustbin should be > 0
[42]	567
[1720]	568	if ((microphys.and..not.doubleq).or.
[358]	569	& (microphys.and..not.water)) then
[1720]	570	print*,'if microphys is used, then doubleq,'
	571	print*,'and water must be used!'
	572	stop
	573	endif
	574	if (microphys.and..not.scavenging) then
	575	print*,''
	576	print*,'----------------WARNING-----------------'
	577	print*,'microphys is used without scavenging !!!'
	578	print*,'----------------WARNING-----------------'
	579	print*,''
	580	endif
	581
	582	if ((scavenging.and..not.microphys).or.
[1617]	583	& (scavenging.and.(dustbin.lt.1)))then
[1720]	584	print*,'if scavenging is used, then microphys'
	585	print*,'must be used!'
	586	stop
	587	endif
[740]	588
[358]	589	! Test of incompatibility:
	590
[42]	591	write(,) "Permanent water caps at poles ?",
	592	& " .true. is RECOMMENDED"
	593	write(,) "(with .true., North cap is a source of water ",
	594	& "and South pole is a cold trap)"
	595	caps=.true. ! default value
	596	call getin("caps",caps)
	597	write(,) " caps = ",caps
	598
[358]	599	! albedo_h2o_ice
[283]	600	write(,) "water ice albedo ?"
	601	albedo_h2o_ice=0.45
	602	call getin("albedo_h2o_ice",albedo_h2o_ice)
	603	write(,) " albedo_h2o_ice = ",albedo_h2o_ice
[358]	604	! inert_h2o_ice
[283]	605	write(,) "water ice thermal inertia ?"
	606	inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
	607	call getin("inert_h2o_ice",inert_h2o_ice)
	608	write(,) " inert_h2o_ice = ",inert_h2o_ice
[358]	609	! frost_albedo_threshold
[283]	610	write(,) "frost thickness threshold for albedo ?"
[285]	611	frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2)
[283]	612	call getin("frost_albedo_threshold",
	613	& frost_albedo_threshold)
	614	write(,) " frost_albedo_threshold = ",
	615	& frost_albedo_threshold
	616
[485]	617	! call Titus crocus line -- DEFAULT IS NONE
	618	write(,) "Titus crocus line ?"
	619	tituscap=.false. ! default value
	620	call getin("tituscap",tituscap)
	621	write(,) "tituscap",tituscap
[633]	622
[485]	623
[42]	624	write(,) "photochemistry: include chemical species"
	625	photochem=.false. ! default value
	626	call getin("photochem",photochem)
	627	write(,) " photochem = ",photochem
	628
	629
[1246]	630	! SCATTERERS
	631	write(,) "how many scatterers?"
	632	naerkind=1 ! default value
	633	call getin("naerkind",naerkind)
	634	write(,)" naerkind = ",naerkind
	635
	636	! Test of incompatibility
	637	c Logical tests for radiatively active water-ice clouds:
	638	IF ( (activice.AND.(.NOT.water)).OR.
	639	& (activice.AND.(naerkind.LT.2)) ) THEN
	640	WRITE(,) 'If activice is TRUE, water has to be set'
	641	WRITE(,) 'to TRUE, and "naerkind" must be at least'
	642	WRITE(,) 'equal to 2.'
	643	CALL ABORT
	644	ELSE IF ( (.NOT.activice).AND.(naerkind.GT.1) ) THEN
	645	WRITE(,) 'naerkind is greater than unity, but'
	646	WRITE(,) 'activice has not been set to .true.'
	647	WRITE(,) 'in callphys.def; this is not logical!'
	648	CALL ABORT
	649	ENDIF
	650
	651	!------------------------------------------
	652	!------------------------------------------
	653	! once naerkind is known allocate arrays
	654	! -- we do it here and not in phys_var_init
	655	! -- because we need to know naerkind first
	656	CALL ini_scatterers(ngrid,nlayer)
	657	!------------------------------------------
	658	!------------------------------------------
	659
	660
	661	c Please name the different scatterers here ----------------
	662	name_iaer(1) = "dust_conrath" !! default choice is good old Conrath profile
	663	IF (doubleq.AND.active) name_iaer(1) = "dust_doubleq" !! two-moment scheme
	664	if (nq.gt.1) then
	665	! trick to avoid problems compiling with 1 tracer
	666	! and picky compilers who know name_iaer(2) is out of bounds
	667	j=2
	668	IF (water.AND.activice) name_iaer(j) = "h2o_ice" !! radiatively-active clouds
	669	IF (submicron.AND.active) name_iaer(j) = "dust_submicron" !! JBM experimental stuff
	670	endif ! of if (nq.gt.1)
	671	c ----------------------------------------------------------
	672
[42]	673	! THERMOSPHERE
	674
	675	write(,) "call thermosphere ?"
	676	callthermos=.false. ! default value
	677	call getin("callthermos",callthermos)
	678	write(,) " callthermos = ",callthermos
	679
	680
	681	write(,) " water included without cycle ",
	682	& "(only if water=.false.)"
	683	thermoswater=.false. ! default value
	684	call getin("thermoswater",thermoswater)
	685	write(,) " thermoswater = ",thermoswater
	686
	687	write(,) "call thermal conduction ?",
	688	& " (only if callthermos=.true.)"
	689	callconduct=.false. ! default value
	690	call getin("callconduct",callconduct)
	691	write(,) " callconduct = ",callconduct
	692
	693	write(,) "call EUV heating ?",
	694	& " (only if callthermos=.true.)"
	695	calleuv=.false. ! default value
	696	call getin("calleuv",calleuv)
	697	write(,) " calleuv = ",calleuv
	698
	699	write(,) "call molecular viscosity ?",
	700	& " (only if callthermos=.true.)"
	701	callmolvis=.false. ! default value
	702	call getin("callmolvis",callmolvis)
	703	write(,) " callmolvis = ",callmolvis
	704
	705	write(,) "call molecular diffusion ?",
	706	& " (only if callthermos=.true.)"
	707	callmoldiff=.false. ! default value
	708	call getin("callmoldiff",callmoldiff)
	709	write(,) " callmoldiff = ",callmoldiff
	710
	711
	712	write(,) "call thermospheric photochemistry ?",
	713	& " (only if callthermos=.true.)"
	714	thermochem=.false. ! default value
	715	call getin("thermochem",thermochem)
	716	write(,) " thermochem = ",thermochem
	717
[705]	718	write(,) "Method to include solar variability"
[1684]	719	write(,) "0-> fixed value of E10.7 (fixed_euv_value); ",
	720	& "1-> daily evolution of E10.7 (for given solvaryear)"
[705]	721	solvarmod=1
	722	call getin("solvarmod",solvarmod)
	723	write(,) " solvarmod = ",solvarmod
	724
[1684]	725	write(,) "Fixed euv (for solvarmod==0) 10.7 value?"
	726	write(,) " (min=80 , ave=140, max=320)"
	727	fixed_euv_value=140 ! default value
	728	call getin("fixed_euv_value",fixed_euv_value)
	729	write(,) " fixed_euv_value = ",fixed_euv_value
[42]	730
[705]	731	write(,) "Solar variability as observed for MY: "
	732	write(,) "Only if solvarmod=1"
	733	solvaryear=24
	734	call getin("solvaryear",solvaryear)
	735	write(,) " solvaryear = ",solvaryear
	736
[552]	737	write(,) "UV heating efficiency:",
	738	& "measured values between 0.19 and 0.23 (Fox et al. 1996)",
	739	& "lower values may be used to compensate low 15 um cooling"
	740	euveff=0.21 !default value
	741	call getin("euveff",euveff)
	742	write(,) " euveff = ", euveff
[42]	743
[1651]	744
[42]	745	if (.not.callthermos) then
	746	if (thermoswater) then
	747	print*,'if thermoswater is set, callthermos must be true'
	748	stop
	749	endif
	750	if (callconduct) then
	751	print*,'if callconduct is set, callthermos must be true'
	752	stop
	753	endif
	754	if (calleuv) then
	755	print*,'if calleuv is set, callthermos must be true'
	756	stop
	757	endif
	758	if (callmolvis) then
	759	print*,'if callmolvis is set, callthermos must be true'
	760	stop
	761	endif
	762	if (callmoldiff) then
	763	print*,'if callmoldiff is set, callthermos must be true'
	764	stop
	765	endif
	766	if (thermochem) then
	767	print*,'if thermochem is set, callthermos must be true'
	768	stop
	769	endif
	770	endif
	771
	772	! Test of incompatibility:
	773	! if photochem is used, then water should be used too
	774
	775	if (photochem.and..not.water) then
	776	print*,'if photochem is used, water should be used too'
	777	stop
	778	endif
	779
	780	! if callthermos is used, then thermoswater should be used too
	781	! (if water not used already)
	782
	783	if (callthermos .and. .not.water) then
	784	if (callthermos .and. .not.thermoswater) then
	785	print*,'if callthermos is used, water or thermoswater
	786	& should be used too'
	787	stop
	788	endif
	789	endif
	790
	791	PRINT*,'--------------------------------------------'
	792	PRINT*
	793	PRINT*
	794	ELSE
	795	write(,)
	796	write(,) 'Cannot read file callphys.def. Is it here ?'
	797	stop
	798	ENDIF
	799
	800	8000 FORMAT(t5,a12,l8)
	801	8001 FORMAT(t5,a12,i8)
	802
	803	PRINT*
[1224]	804	PRINT*,'conf_phys: daysec',daysec
[42]	805	PRINT*
[1224]	806	PRINT*,'conf_phys: The radiative transfer is computed:'
[42]	807	PRINT*,' each ',iradia,' physical time-step'
	808	PRINT,' or each ',iradiadtphys,' seconds'
	809	PRINT*
	810	! --------------------------------------------------------------
	811	! Managing the Longwave radiative transfer
	812	! --------------------------------------------------------------
	813
	814	! In most cases, the run just use the following values :
	815	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	816	callemis=.true.
	817	! ilwd=10*int(daysec/dtphys) ! bug before 22/10/01
	818	ilwd=1
	819	ilwn=1 !2
	820	ilwb=1 !2
	821	linear=.true.
	822	ncouche=3
	823	alphan=0.4
	824	semi=0
	825
	826	! BUT people working hard on the LW may want to read them in 'radia.def'
	827	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	828
	829	OPEN(99,file='radia.def',status='old',form='formatted'
	830	. ,iostat=ierr)
	831	IF(ierr.EQ.0) THEN
[1224]	832	write(,) 'conf_phys: Reading radia.def !!!'
[42]	833	READ(99,fmt='(a)') ch1
	834	READ(99,*) callemis
	835	WRITE(*,8000) ch1,callemis
	836
	837	READ(99,fmt='(a)') ch1
	838	READ(99,*) iradia
	839	WRITE(*,8001) ch1,iradia
	840
	841	READ(99,fmt='(a)') ch1
	842	READ(99,*) ilwd
	843	WRITE(*,8001) ch1,ilwd
	844
	845	READ(99,fmt='(a)') ch1
	846	READ(99,*) ilwn
	847	WRITE(*,8001) ch1,ilwn
	848
	849	READ(99,fmt='(a)') ch1
	850	READ(99,*) linear
	851	WRITE(*,8000) ch1,linear
	852
	853	READ(99,fmt='(a)') ch1
	854	READ(99,*) ncouche
	855	WRITE(*,8001) ch1,ncouche
	856
	857	READ(99,fmt='(a)') ch1
	858	READ(99,*) alphan
	859	WRITE(,) ch1,alphan
	860
	861	READ(99,fmt='(a)') ch1
	862	READ(99,*) ilwb
	863	WRITE(*,8001) ch1,ilwb
	864
	865
	866	READ(99,fmt='(a)') ch1
	867	READ(99,'(l1)') callg2d
	868	WRITE(*,8000) ch1,callg2d
	869
	870	READ(99,fmt='(a)') ch1
	871	READ(99,*) semi
	872	WRITE(,) ch1,semi
	873	end if
	874	CLOSE(99)
	875
	876	END

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

Download in other formats:

Original Format