Context Navigation

conf_phys.F @ 1918

Last change on this file since 1918 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

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	use surfdat_h, only: albedo_h2o_ice, inert_h2o_ice,
40	& frost_albedo_threshold
41	use time_phylmdz_mod, only: ecritphy,day_step,iphysiq,ecritstart,
42	& daysec,dtphys
43	use dimradmars_mod, only: naerkind, name_iaer,
44	& ini_scatterers,tauvis
45	use datafile_mod, only: datadir
46
47	IMPLICIT NONE
48	include "callkeys.h"
49	include "microphys.h"
50
51	INTEGER,INTENT(IN) :: ngrid,nlayer,nq
52	INTEGER ierr,j
53
54	CHARACTER ch1*12
55	#ifndef MESOSCALE
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
67	#endif
68
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*,'--------------------------------------------'
82	PRINT*,' conf_phys: Parameters for the physics (callphys.def)'
83	PRINT*,'--------------------------------------------'
84
85	write(,) "Directory where external input files are:"
86	! default path is set in datafile_mod
87	call getin("datadir",datadir)
88	write(,) " datadir = ",trim(datadir)
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 ?"
114	#ifdef MESOSCALE
115	callstats=.false. ! default value
116	#else
117	callstats=.true. ! default value
118	#endif
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,",
130	& "=6 cold (low dust) scenario; =7 warm (high dust) scenario ",
131	& "=24,25 ... 30 :Mars Year 24, ... or 30 from TES assimilation"
132	iaervar=3 ! default value
133	call getin("iaervar",iaervar)
134	write(,) " iaervar = ",iaervar
135
136	write(,) "Reference (visible) dust opacity at 610 Pa ",
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
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
160	write(,) "call radiative transfer ?"
161	callrad=.true. ! default value
162	call getin("callrad",callrad)
163	write(,) " callrad = ",callrad
164
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
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
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
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
195
196	write(,) "New NIR NLTE correction ?",
197	$ "0-> old model (no correction)",
198	$ "1-> new correction",
199	$ "(matters only if callnirco2=T)"
200	#ifdef MESOSCALE
201	nircorr=0 !default value. this is OK below 60 km.
202	#else
203	nircorr=0 !default value
204	#endif
205	call getin("nircorr",nircorr)
206	write(,) " nircorr = ",nircorr
207
208	write(,) "call turbulent vertical diffusion ?"
209	calldifv=.true. ! default value
210	call getin("calldifv",calldifv)
211	write(,) " calldifv = ",calldifv
212
213	write(,) "call thermals ?"
214	calltherm=.false. ! default value
215	call getin("calltherm",calltherm)
216	write(,) " calltherm = ",calltherm
217
218	write(,) "call convective adjustment ?"
219	calladj=.true. ! default value
220	call getin("calladj",calladj)
221	write(,) " calladj = ",calladj
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'
228	endif
229
230	write(,) "used latest version of yamada scheme?"
231	callyamada4=.true. ! default value
232	call getin("callyamada4",callyamada4)
233	write(,) " callyamada4 = ",callyamada4
234
235	if (calltherm .and. .not.callyamada4) then
236	print*,'!!!! WARNING WARNING WARNING !!!!'
237	print*,'if calltherm=T we strongly advise that '
238	print*,'you set the flag callyamada4 to T '
239	print*,'!!!! WARNING WARNING WARNING !!!!'
240	endif
241
242	write(,) "call Richardson-based surface layer ?"
243	callrichsl=.false. ! default value
244	call getin("callrichsl",callrichsl)
245	write(,) " callrichsl = ",callrichsl
246
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
254	if (calladj .and. callrichsl .and. (.not. calltherm)) then
255	print*,'You should not be calling the convective adjustment
256	& scheme with the Richardson surface-layer and without the thermals
257	&. This approach is not
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
264
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
304	! dustbin
305	write(,)"Transported dust ? (if >0, use 'dustbin' dust bins)"
306	dustbin=0 ! default value
307	call getin("dustbin",dustbin)
308	write(,)" dustbin = ",dustbin
309	! active
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
322	! doubleq
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
328	! submicron
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
340	if (doubleq.and.submicron.and.(nq.LT.3)) then
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
346
347	! lifting
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
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
371	#ifndef MESOSCALE
372	! this test is valid in GCM case
373	! ... not in mesoscale case, for we want to activate mesoscale lifting
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
379	#endif
380
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
407	! callddevil
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
420	! sedimentation
421	write(,) "Gravitationnal sedimentation ?"
422	sedimentation=.true. ! default value
423	call getin("sedimentation",sedimentation)
424	write(,) " sedimentation = ",sedimentation
425	! activice
426	write(,) "Radiatively active transported atmospheric ",
427	& "water ice ?"
428	activice=.false. ! default value
429	call getin("activice",activice)
430	write(,) " activice = ",activice
431	! water
432	write(,) "Compute water cycle ?"
433	water=.false. ! default value
434	call getin("water",water)
435	write(,) " water = ",water
436
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
448	!CO2 clouds scheme?
449	write(,) "Compute CO2 clouds (implies microphysical scheme)?"
450	co2clouds=.false. ! default value
451	call getin("co2clouds",co2clouds)
452	write(,) " co2clouds = ",co2clouds
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
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
478
479
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
486	! Test of incompatibility:
487
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
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
509
510	! water ice clouds effective variance distribution for sedimentaion
511	write(,) "Sed effective variance for water ice clouds ?"
512	nuice_sed=0.45
513	call getin("nuice_sed",nuice_sed)
514	write(,) "water_param nueff Sedimentation:", nuice_sed
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
532	! ccn factor if no scavenging
533	write(,) "water param CCN reduc. factor ?"
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
540	! microphys
541	write(,)"Microphysical scheme for water-ice clouds?"
542	microphys=.false. ! default value
543	call getin("microphys",microphys)
544	write(,)" microphys = ",microphys
545
546	! supersat
547	write(,)"Allow super-saturation of water vapor?"
548	supersat=.true. ! default value
549	call getin("supersat",supersat)
550	write(,)"supersat = ",supersat
551
552	! microphysical parameter contact
553	write(,) "water contact parameter ?"
554	mteta = 0.95
555	call getin("mteta",mteta)
556	write(,) "mteta = ", mteta
557
558	! scavenging
559	write(,)"Dust scavenging by H2O/CO2 snowfall ?"
560	scavenging=.false. ! default value
561	call getin("scavenging",scavenging)
562	write(,)" scavenging = ",scavenging
563
564
565	! Test of incompatibility:
566	! if scavenging is used, then dustbin should be > 0
567
568	if ((microphys.and..not.doubleq).or.
569	& (microphys.and..not.water)) then
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.
583	& (scavenging.and.(dustbin.lt.1)))then
584	print*,'if scavenging is used, then microphys'
585	print*,'must be used!'
586	stop
587	endif
588
589	! Test of incompatibility:
590
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
599	! albedo_h2o_ice
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
604	! inert_h2o_ice
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
609	! frost_albedo_threshold
610	write(,) "frost thickness threshold for albedo ?"
611	frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2)
612	call getin("frost_albedo_threshold",
613	& frost_albedo_threshold)
614	write(,) " frost_albedo_threshold = ",
615	& frost_albedo_threshold
616
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
622
623
624	write(,) "photochemistry: include chemical species"
625	photochem=.false. ! default value
626	call getin("photochem",photochem)
627	write(,) " photochem = ",photochem
628
629
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
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
718	write(,) "Method to include solar variability"
719	write(,) "0-> fixed value of E10.7 (fixed_euv_value); ",
720	& "1-> daily evolution of E10.7 (for given solvaryear)"
721	solvarmod=1
722	call getin("solvarmod",solvarmod)
723	write(,) " solvarmod = ",solvarmod
724
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
730
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
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
743
744
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*
804	PRINT*,'conf_phys: daysec',daysec
805	PRINT*
806	PRINT*,'conf_phys: The radiative transfer is computed:'
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
832	write(,) 'conf_phys: Reading radia.def !!!'
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