Context Navigation

conf_phys.F @ 1264

Last change on this file since 1264 was 1264, checked in by aslmd, 11 years ago
LMDZ.MARS. Two modifications related to outputs. 1. corrected a bug for thermospheric outputs (writediagfi under a callstats flag) 2. (only for mesoscale) made adaptations to allow for freedust runs with MESOSCALE.
File size: 26.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
[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 "dimensions.h"
	48	#include "dimphys.h"
[1047]	49	!#include "comsaison.h"
	50	!#include "comdiurn.h"
	51	!#include "comgeomfi.h"
[42]	52	#include "callkeys.h"
[1047]	53	!#include "surfdat.h"
	54	!#include "dimradmars.h"
[42]	55	#include "datafile.h"
[1047]	56	!#include "slope.h"
[358]	57	#include "microphys.h"
[1036]	58	!#include "tracer.h"
[1112]	59
[1246]	60	INTEGER,INTENT(IN) :: ngrid,nlayer,nq
	61	INTEGER ig,ierr,j
[42]	62
	63	CHARACTER ch1*12
	64	CHARACTER ch80*80
	65
[1130]	66	! read in 'ecritphy' (frequency of calls to physics, in dynamical steps)
	67	! (also done in dyn3d/defrun_new but not in LMDZ.COMMON)
	68	call getin("ecritphy",ecritphy)
	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:"
[426]	87	datafile="/u/forget/WWW/datagcm/datafile"
[226]	88	call getin("datadir",datafile) ! default path
[42]	89	write(,) " datafile = ",trim(datafile)
	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
	283
	284	write(,)"rad.transfer is computed every iradia",
	285	& " physical timestep"
	286	iradia=1 ! default value
	287	call getin("iradia",iradia)
	288	write(,)" iradia = ",iradia
	289
	290
	291	write(,)"Output of the exchange coefficient mattrix ?",
	292	& "(for diagnostics only)"
	293	callg2d=.false. ! default value
	294	call getin("callg2d",callg2d)
	295	write(,)" callg2d = ",callg2d
	296
	297	write(,)"Rayleigh scattering : (should be .false. for now)"
	298	rayleigh=.false.
	299	call getin("rayleigh",rayleigh)
	300	write(,)" rayleigh = ",rayleigh
	301
	302
	303	! TRACERS:
	304
[358]	305	! dustbin
[42]	306	write(,)"Transported dust ? (if >0, use 'dustbin' dust bins)"
	307	dustbin=0 ! default value
	308	call getin("dustbin",dustbin)
	309	write(,)" dustbin = ",dustbin
[358]	310	! active
[42]	311	write(,)"Radiatively active dust ? (matters if dustbin>0)"
	312	active=.false. ! default value
	313	call getin("active",active)
	314	write(,)" active = ",active
	315
	316	! Test of incompatibility:
	317	! if active is used, then dustbin should be > 0
	318
	319	if (active.and.(dustbin.lt.1)) then
	320	print*,'if active is used, then dustbin should > 0'
	321	stop
	322	endif
[358]	323	! doubleq
[42]	324	write(,)"use mass and number mixing ratios to predict",
	325	& " dust size ?"
	326	doubleq=.false. ! default value
	327	call getin("doubleq",doubleq)
	328	write(,)" doubleq = ",doubleq
[358]	329	! submicron
[42]	330	submicron=.false. ! default value
	331	call getin("submicron",submicron)
	332	write(,)" submicron = ",submicron
	333
	334	! Test of incompatibility:
	335	! if doubleq is used, then dustbin should be 2
	336
	337	if (doubleq.and.(dustbin.ne.2)) then
	338	print*,'if doubleq is used, then dustbin should be 2'
	339	stop
	340	endif
[1036]	341	if (doubleq.and.submicron.and.(nq.LT.3)) then
[42]	342	print*,'If doubleq is used with a submicron tracer,'
	343	print*,' then the number of tracers has to be'
	344	print*,' larger than 3.'
	345	stop
	346	endif
[1088]	347
[358]	348	! lifting
[42]	349	write(,)"dust lifted by GCM surface winds ?"
	350	lifting=.false. ! default value
	351	call getin("lifting",lifting)
	352	write(,)" lifting = ",lifting
	353
	354	! Test of incompatibility:
	355	! if lifting is used, then dustbin should be > 0
	356
	357	if (lifting.and.(dustbin.lt.1)) then
	358	print*,'if lifting is used, then dustbin should > 0'
	359	stop
	360	endif
[1088]	361
	362	! free evolving dust
	363	! freedust=true just says that there is no lifting and no dust opacity scaling.
	364	write(,)"dust lifted by GCM surface winds ?"
	365	freedust=.false. ! default value
	366	call getin("freedust",freedust)
	367	write(,)" freedust = ",freedust
	368	if (freedust.and..not.doubleq) then
	369	print*,'freedust should be used with doubleq !'
	370	stop
	371	endif
[1264]	372	#ifndef MESOSCALE
	373	! this test is valid in GCM case
	374	! ... not in mesoscale case, for we want to activate mesoscale lifting
[1088]	375	if (freedust.and.lifting) then
	376	print*,'if freedust is used, then lifting should not be used'
	377	print*,'lifting forced to false !!'
	378	lifting=.false.
	379	endif
[1264]	380	#endif
[1088]	381
[358]	382	! callddevil
[42]	383	write(,)" dust lifted by dust devils ?"
	384	callddevil=.false. !default value
	385	call getin("callddevil",callddevil)
	386	write(,)" callddevil = ",callddevil
	387
	388	! Test of incompatibility:
	389	! if dustdevil is used, then dustbin should be > 0
	390
	391	if (callddevil.and.(dustbin.lt.1)) then
	392	print*,'if dustdevil is used, then dustbin should > 0'
	393	stop
	394	endif
[358]	395	! sedimentation
[42]	396	write(,) "Gravitationnal sedimentation ?"
	397	sedimentation=.true. ! default value
	398	call getin("sedimentation",sedimentation)
	399	write(,) " sedimentation = ",sedimentation
[358]	400	! activice
[42]	401	write(,) "Radiatively active transported atmospheric ",
	402	& "water ice ?"
	403	activice=.false. ! default value
	404	call getin("activice",activice)
	405	write(,) " activice = ",activice
[358]	406	! water
[42]	407	write(,) "Compute water cycle ?"
	408	water=.false. ! default value
	409	call getin("water",water)
	410	write(,) " water = ",water
	411
[833]	412	! thermal inertia feedback
	413	write(,) "Activate the thermal inertia feedback ?"
	414	tifeedback=.false. ! default value
	415	call getin("tifeedback",tifeedback)
	416	write(,) " tifeedback = ",tifeedback
	417
[42]	418	! Test of incompatibility:
	419
[833]	420	if (tifeedback.and..not.water) then
	421	print*,'if tifeedback is used,'
	422	print*,'water should be used too'
	423	stop
	424	endif
	425
	426	if (tifeedback.and..not.callsoil) then
	427	print*,'if tifeedback is used,'
	428	print*,'callsoil should be used too'
	429	stop
	430	endif
	431
[42]	432	if (activice.and..not.water) then
	433	print*,'if activice is used, water should be used too'
	434	stop
	435	endif
	436
	437	if (water.and..not.tracer) then
	438	print*,'if water is used, tracer should be used too'
	439	stop
	440	endif
[420]	441
[455]	442	! water ice clouds effective variance distribution for sedimentaion
	443	write(,) "effective variance for water ice clouds ?"
	444	nuice_sed=0.45
	445	call getin("nuice_sed",nuice_sed)
	446	write(,) "water_param nueff Sedimentation:", nuice_sed
	447
[420]	448	! ccn factor if no scavenging
	449	write(,) "water param CCN reduc. factor ?", ccn_factor
	450	ccn_factor = 4.5
	451	call getin("ccn_factor",ccn_factor)
	452	write(,)" ccn_factor = ",ccn_factor
	453	write(,)"Careful: only used when microphys=F, otherwise"
	454	write(,)"the contact parameter is used instead;"
	455
[358]	456	! microphys
	457	write(,)"Microphysical scheme for water-ice clouds?"
	458	microphys=.false. ! default value
	459	call getin("microphys",microphys)
	460	write(,)" microphys = ",microphys
[42]	461
[358]	462	! microphysical parameter contact
	463	write(,) "water contact parameter ?"
	464	mteta = 0.95
	465	call getin("mteta",mteta)
	466	write(,) "mteta = ", mteta
	467
	468	! scavenging
	469	write(,)"Dust scavenging by H2O/CO2 snowfall ?"
	470	scavenging=.false. ! default value
	471	call getin("scavenging",scavenging)
	472	write(,)" scavenging = ",scavenging
	473
	474
[42]	475	! Test of incompatibility:
[358]	476	! if scavenging is used, then dustbin should be > 0
[42]	477
[358]	478	if ((microphys.and..not.doubleq).or.
	479	& (microphys.and..not.water)) then
[740]	480	print*,'if microphys is used, then doubleq,'
	481	print*,'and water must be used!'
	482	stop
[358]	483	endif
[740]	484	if (microphys.and..not.scavenging) then
	485	print*,''
	486	print*,'----------------WARNING-----------------'
	487	print*,'microphys is used without scavenging !!!'
	488	print*,'----------------WARNING-----------------'
	489	print*,''
[358]	490	endif
	491
[740]	492	if ((scavenging.and..not.microphys).or.
	493	& (scavenging.and.(dustbin.lt.1))) then
	494	print*,'if scavenging is used, then microphys'
	495	print*,'must be used!'
	496	stop
	497	endif
	498
[358]	499	! Test of incompatibility:
	500
[42]	501	write(,) "Permanent water caps at poles ?",
	502	& " .true. is RECOMMENDED"
	503	write(,) "(with .true., North cap is a source of water ",
	504	& "and South pole is a cold trap)"
	505	caps=.true. ! default value
	506	call getin("caps",caps)
	507	write(,) " caps = ",caps
	508
[358]	509	! albedo_h2o_ice
[283]	510	write(,) "water ice albedo ?"
	511	albedo_h2o_ice=0.45
	512	call getin("albedo_h2o_ice",albedo_h2o_ice)
	513	write(,) " albedo_h2o_ice = ",albedo_h2o_ice
[358]	514	! inert_h2o_ice
[283]	515	write(,) "water ice thermal inertia ?"
	516	inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
	517	call getin("inert_h2o_ice",inert_h2o_ice)
	518	write(,) " inert_h2o_ice = ",inert_h2o_ice
[358]	519	! frost_albedo_threshold
[283]	520	write(,) "frost thickness threshold for albedo ?"
[285]	521	frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2)
[283]	522	call getin("frost_albedo_threshold",
	523	& frost_albedo_threshold)
	524	write(,) " frost_albedo_threshold = ",
	525	& frost_albedo_threshold
	526
[485]	527	! call Titus crocus line -- DEFAULT IS NONE
	528	write(,) "Titus crocus line ?"
	529	tituscap=.false. ! default value
	530	call getin("tituscap",tituscap)
	531	write(,) "tituscap",tituscap
[633]	532
[485]	533
[42]	534	write(,) "photochemistry: include chemical species"
	535	photochem=.false. ! default value
	536	call getin("photochem",photochem)
	537	write(,) " photochem = ",photochem
	538
	539
[1246]	540	! SCATTERERS
	541	write(,) "how many scatterers?"
	542	naerkind=1 ! default value
	543	call getin("naerkind",naerkind)
	544	write(,)" naerkind = ",naerkind
	545
	546	! Test of incompatibility
	547	c Logical tests for radiatively active water-ice clouds:
	548	IF ( (activice.AND.(.NOT.water)).OR.
	549	& (activice.AND.(naerkind.LT.2)) ) THEN
	550	WRITE(,) 'If activice is TRUE, water has to be set'
	551	WRITE(,) 'to TRUE, and "naerkind" must be at least'
	552	WRITE(,) 'equal to 2.'
	553	CALL ABORT
	554	ELSE IF ( (.NOT.activice).AND.(naerkind.GT.1) ) THEN
	555	WRITE(,) 'naerkind is greater than unity, but'
	556	WRITE(,) 'activice has not been set to .true.'
	557	WRITE(,) 'in callphys.def; this is not logical!'
	558	CALL ABORT
	559	ENDIF
	560
	561	!------------------------------------------
	562	!------------------------------------------
	563	! once naerkind is known allocate arrays
	564	! -- we do it here and not in phys_var_init
	565	! -- because we need to know naerkind first
	566	CALL ini_scatterers(ngrid,nlayer)
	567	!------------------------------------------
	568	!------------------------------------------
	569
	570
	571	c Please name the different scatterers here ----------------
	572	name_iaer(1) = "dust_conrath" !! default choice is good old Conrath profile
	573	IF (doubleq.AND.active) name_iaer(1) = "dust_doubleq" !! two-moment scheme
	574	if (nq.gt.1) then
	575	! trick to avoid problems compiling with 1 tracer
	576	! and picky compilers who know name_iaer(2) is out of bounds
	577	j=2
	578	IF (water.AND.activice) name_iaer(j) = "h2o_ice" !! radiatively-active clouds
	579	IF (submicron.AND.active) name_iaer(j) = "dust_submicron" !! JBM experimental stuff
	580	endif ! of if (nq.gt.1)
	581	c ----------------------------------------------------------
	582
[42]	583	! THERMOSPHERE
	584
	585	write(,) "call thermosphere ?"
	586	callthermos=.false. ! default value
	587	call getin("callthermos",callthermos)
	588	write(,) " callthermos = ",callthermos
	589
	590
	591	write(,) " water included without cycle ",
	592	& "(only if water=.false.)"
	593	thermoswater=.false. ! default value
	594	call getin("thermoswater",thermoswater)
	595	write(,) " thermoswater = ",thermoswater
	596
	597	write(,) "call thermal conduction ?",
	598	& " (only if callthermos=.true.)"
	599	callconduct=.false. ! default value
	600	call getin("callconduct",callconduct)
	601	write(,) " callconduct = ",callconduct
	602
	603	write(,) "call EUV heating ?",
	604	& " (only if callthermos=.true.)"
	605	calleuv=.false. ! default value
	606	call getin("calleuv",calleuv)
	607	write(,) " calleuv = ",calleuv
	608
	609	write(,) "call molecular viscosity ?",
	610	& " (only if callthermos=.true.)"
	611	callmolvis=.false. ! default value
	612	call getin("callmolvis",callmolvis)
	613	write(,) " callmolvis = ",callmolvis
	614
	615	write(,) "call molecular diffusion ?",
	616	& " (only if callthermos=.true.)"
	617	callmoldiff=.false. ! default value
	618	call getin("callmoldiff",callmoldiff)
	619	write(,) " callmoldiff = ",callmoldiff
	620
	621
	622	write(,) "call thermospheric photochemistry ?",
	623	& " (only if callthermos=.true.)"
	624	thermochem=.false. ! default value
	625	call getin("thermochem",thermochem)
	626	write(,) " thermochem = ",thermochem
	627
[705]	628	write(,) "Method to include solar variability"
	629	write(,) "0-> old method (using solarcondate); ",
	630	& "1-> variability wit E10.7"
	631	solvarmod=1
	632	call getin("solvarmod",solvarmod)
	633	write(,) " solvarmod = ",solvarmod
	634
[42]	635	write(,) "date for solar flux calculation:",
[705]	636	& " (1985 < date < 2002)",
	637	$ " (Only used if solvarmod=0)"
[42]	638	write(,) "(Solar min=1996.4 ave=1993.4 max=1990.6)"
	639	solarcondate=1993.4 ! default value
	640	call getin("solarcondate",solarcondate)
	641	write(,) " solarcondate = ",solarcondate
	642
[705]	643	write(,) "Solar variability as observed for MY: "
	644	write(,) "Only if solvarmod=1"
	645	solvaryear=24
	646	call getin("solvaryear",solvaryear)
	647	write(,) " solvaryear = ",solvaryear
	648
[552]	649	write(,) "UV heating efficiency:",
	650	& "measured values between 0.19 and 0.23 (Fox et al. 1996)",
	651	& "lower values may be used to compensate low 15 um cooling"
	652	euveff=0.21 !default value
	653	call getin("euveff",euveff)
	654	write(,) " euveff = ", euveff
[42]	655
	656	if (.not.callthermos) then
	657	if (thermoswater) then
	658	print*,'if thermoswater is set, callthermos must be true'
	659	stop
	660	endif
	661	if (callconduct) then
	662	print*,'if callconduct is set, callthermos must be true'
	663	stop
	664	endif
	665	if (calleuv) then
	666	print*,'if calleuv is set, callthermos must be true'
	667	stop
	668	endif
	669	if (callmolvis) then
	670	print*,'if callmolvis is set, callthermos must be true'
	671	stop
	672	endif
	673	if (callmoldiff) then
	674	print*,'if callmoldiff is set, callthermos must be true'
	675	stop
	676	endif
	677	if (thermochem) then
	678	print*,'if thermochem is set, callthermos must be true'
	679	stop
	680	endif
	681	endif
	682
	683	! Test of incompatibility:
	684	! if photochem is used, then water should be used too
	685
	686	if (photochem.and..not.water) then
	687	print*,'if photochem is used, water should be used too'
	688	stop
	689	endif
	690
	691	! if callthermos is used, then thermoswater should be used too
	692	! (if water not used already)
	693
	694	if (callthermos .and. .not.water) then
	695	if (callthermos .and. .not.thermoswater) then
	696	print*,'if callthermos is used, water or thermoswater
	697	& should be used too'
	698	stop
	699	endif
	700	endif
	701
	702	PRINT*,'--------------------------------------------'
	703	PRINT*
	704	PRINT*
	705	ELSE
	706	write(,)
	707	write(,) 'Cannot read file callphys.def. Is it here ?'
	708	stop
	709	ENDIF
	710
	711	8000 FORMAT(t5,a12,l8)
	712	8001 FORMAT(t5,a12,i8)
	713
	714	PRINT*
[1224]	715	PRINT*,'conf_phys: daysec',daysec
[42]	716	PRINT*
[1224]	717	PRINT*,'conf_phys: The radiative transfer is computed:'
[42]	718	PRINT*,' each ',iradia,' physical time-step'
	719	PRINT,' or each ',iradiadtphys,' seconds'
	720	PRINT*
	721	! --------------------------------------------------------------
	722	! Managing the Longwave radiative transfer
	723	! --------------------------------------------------------------
	724
	725	! In most cases, the run just use the following values :
	726	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	727	callemis=.true.
	728	! ilwd=10*int(daysec/dtphys) ! bug before 22/10/01
	729	ilwd=1
	730	ilwn=1 !2
	731	ilwb=1 !2
	732	linear=.true.
	733	ncouche=3
	734	alphan=0.4
	735	semi=0
	736
	737	! BUT people working hard on the LW may want to read them in 'radia.def'
	738	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	739
	740	OPEN(99,file='radia.def',status='old',form='formatted'
	741	. ,iostat=ierr)
	742	IF(ierr.EQ.0) THEN
[1224]	743	write(,) 'conf_phys: Reading radia.def !!!'
[42]	744	READ(99,fmt='(a)') ch1
	745	READ(99,*) callemis
	746	WRITE(*,8000) ch1,callemis
	747
	748	READ(99,fmt='(a)') ch1
	749	READ(99,*) iradia
	750	WRITE(*,8001) ch1,iradia
	751
	752	READ(99,fmt='(a)') ch1
	753	READ(99,*) ilwd
	754	WRITE(*,8001) ch1,ilwd
	755
	756	READ(99,fmt='(a)') ch1
	757	READ(99,*) ilwn
	758	WRITE(*,8001) ch1,ilwn
	759
	760	READ(99,fmt='(a)') ch1
	761	READ(99,*) linear
	762	WRITE(*,8000) ch1,linear
	763
	764	READ(99,fmt='(a)') ch1
	765	READ(99,*) ncouche
	766	WRITE(*,8001) ch1,ncouche
	767
	768	READ(99,fmt='(a)') ch1
	769	READ(99,*) alphan
	770	WRITE(,) ch1,alphan
	771
	772	READ(99,fmt='(a)') ch1
	773	READ(99,*) ilwb
	774	WRITE(*,8001) ch1,ilwb
	775
	776
	777	READ(99,fmt='(a)') ch1
	778	READ(99,'(l1)') callg2d
	779	WRITE(*,8000) ch1,callg2d
	780
	781	READ(99,fmt='(a)') ch1
	782	READ(99,*) semi
	783	WRITE(,) ch1,semi
	784	end if
	785	CLOSE(99)
	786
	787	END

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

Download in other formats:

Original Format