Context Navigation

conf_phys.F @ 1353

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

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

Download in other formats:

Original Format