Context Navigation

conf_phys.F @ 1655

Last change on this file since 1655 was 1655, checked in by aslmd, 8 years ago
corrected bugs from rev 1651 that prevented the mesoscale model from compiling. getin does not take double precision arguments. removed meteoriticflux_mod and included the variables in tracer_mod.
File size: 30.3 KB

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

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

Download in other formats:

Original Format