Context Navigation

conf_phys.F @ 1776

Last change on this file since 1776 was 1720, checked in by jaudouard, 7 years ago
Update on CO2 ice clouds scheme
File size: 30.6 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 planete_h
44	use dimradmars_mod, only: naerkind, name_iaer,
45	& ini_scatterers,tauvis
46
47	IMPLICIT NONE
48	#include "callkeys.h"
49	#include "datafile.h"
50	#include "microphys.h"
51
52	INTEGER,INTENT(IN) :: ngrid,nlayer,nq
53	INTEGER ierr,j
54
55	CHARACTER ch1*12
56	#ifndef MESOSCALE
57	! read in some parameters from "run.def" for physics,
58	! or shared between dynamics and physics.
59	ecritphy=240 ! default value
60	call getin("ecritphy",ecritphy) ! frequency of outputs in physics,
61	! in dynamical steps
62	day_step=960 ! default value
63	call getin("day_step",day_step) ! number of dynamical steps per day
64	iphysiq=20 ! default value
65	call getin("iphysiq",iphysiq) ! call physics every iphysiq dyn step
66	ecritstart=0 ! default value
67	call getin("ecritstart",ecritstart) ! write a restart every ecristart steps
68	#endif
69
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*,'--------------------------------------------'
83	PRINT*,' conf_phys: Parameters for the physics (callphys.def)'
84	PRINT*,'--------------------------------------------'
85
86	write(,) "Directory where external input files are:"
87	datafile="/u/lmdz/WWW/planets/mars/datadir"
88	call getin("datadir",datafile) ! default path
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 ?"
115	#ifdef MESOSCALE
116	callstats=.false. ! default value
117	#else
118	callstats=.true. ! default value
119	#endif
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,",
131	& "=6 cold (low dust) scenario; =7 warm (high dust) scenario ",
132	& "=24,25 ... 30 :Mars Year 24, ... or 30 from TES assimilation"
133	iaervar=3 ! default value
134	call getin("iaervar",iaervar)
135	write(,) " iaervar = ",iaervar
136
137	write(,) "Reference (visible) dust opacity at 610 Pa ",
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
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
161	write(,) "call radiative transfer ?"
162	callrad=.true. ! default value
163	call getin("callrad",callrad)
164	write(,) " callrad = ",callrad
165
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
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
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
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
196
197	write(,) "New NIR NLTE correction ?",
198	$ "0-> old model (no correction)",
199	$ "1-> new correction",
200	$ "(matters only if callnirco2=T)"
201	#ifdef MESOSCALE
202	nircorr=0 !default value. this is OK below 60 km.
203	#else
204	nircorr=0 !default value
205	#endif
206	call getin("nircorr",nircorr)
207	write(,) " nircorr = ",nircorr
208
209	write(,) "call turbulent vertical diffusion ?"
210	calldifv=.true. ! default value
211	call getin("calldifv",calldifv)
212	write(,) " calldifv = ",calldifv
213
214	write(,) "call thermals ?"
215	calltherm=.false. ! default value
216	call getin("calltherm",calltherm)
217	write(,) " calltherm = ",calltherm
218
219	write(,) "call convective adjustment ?"
220	calladj=.true. ! default value
221	call getin("calladj",calladj)
222	write(,) " calladj = ",calladj
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'
229	endif
230
231	write(,) "used latest version of yamada scheme?"
232	callyamada4=.true. ! default value
233	call getin("callyamada4",callyamada4)
234	write(,) " callyamada4 = ",callyamada4
235
236	if (calltherm .and. .not.callyamada4) then
237	print*,'!!!! WARNING WARNING WARNING !!!!'
238	print*,'if calltherm=T we strongly advise that '
239	print*,'you set the flag callyamada4 to T '
240	print*,'!!!! WARNING WARNING WARNING !!!!'
241	endif
242
243	write(,) "call Richardson-based surface layer ?"
244	callrichsl=.false. ! default value
245	call getin("callrichsl",callrichsl)
246	write(,) " callrichsl = ",callrichsl
247
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
255	if (calladj .and. callrichsl .and. (.not. calltherm)) then
256	print*,'You should not be calling the convective adjustment
257	& scheme with the Richardson surface-layer and without the thermals
258	&. This approach is not
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
265
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
305	! dustbin
306	write(,)"Transported dust ? (if >0, use 'dustbin' dust bins)"
307	dustbin=0 ! default value
308	call getin("dustbin",dustbin)
309	write(,)" dustbin = ",dustbin
310	! active
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
323	! doubleq
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
329	! submicron
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
341	if (doubleq.and.submicron.and.(nq.LT.3)) then
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
347
348	! lifting
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
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
372	#ifndef MESOSCALE
373	! this test is valid in GCM case
374	! ... not in mesoscale case, for we want to activate mesoscale lifting
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
380	#endif
381
382	! Dust IR opacity
383	write(,)" Wavelength for infrared opacity of dust ?"
384	write(,)" Choices are:"
385	write(,)" tes --- > 9.3 microns [default]"
386	write(,)" mcs --- > 21.6 microns"
387	!
388	! WARNING WARNING WARNING WARNING WARNING WARNING
389	!
390	! BEFORE ADDING A NEW VALUE, BE SURE THAT THE
391	! CORRESPONDING WAVELENGTH IS IN THE LOOKUP TABLE,
392	! OR AT LEAST NO TO FAR, TO AVOID FALLACIOUS INTERPOLATIONS.
393	!
394	dustiropacity="tes" !default value - is expected to shift to mcs one day
395	call getin("dustiropacity",dustiropacity)
396	write(,)" dustiropacity = ",trim(dustiropacity)
397	select case (trim(dustiropacity))
398	case ("tes")
399	dustrefir = 9.3E-6
400	case ("mcs")
401	dustrefir = 21.6E-6
402	case default
403	write(,) trim(dustiropacity),
404	& " is not a valid option for dustiropacity"
405	stop
406	end select
407
408	! callddevil
409	write(,)" dust lifted by dust devils ?"
410	callddevil=.false. !default value
411	call getin("callddevil",callddevil)
412	write(,)" callddevil = ",callddevil
413
414	! Test of incompatibility:
415	! if dustdevil is used, then dustbin should be > 0
416
417	if (callddevil.and.(dustbin.lt.1)) then
418	print*,'if dustdevil is used, then dustbin should > 0'
419	stop
420	endif
421	! sedimentation
422	write(,) "Gravitationnal sedimentation ?"
423	sedimentation=.true. ! default value
424	call getin("sedimentation",sedimentation)
425	write(,) " sedimentation = ",sedimentation
426	! activice
427	write(,) "Radiatively active transported atmospheric ",
428	& "water ice ?"
429	activice=.false. ! default value
430	call getin("activice",activice)
431	write(,) " activice = ",activice
432	! water
433	write(,) "Compute water cycle ?"
434	water=.false. ! default value
435	call getin("water",water)
436	write(,) " water = ",water
437
438	! sub-grid cloud fraction: fixed
439	write(,) "Fixed cloud fraction?"
440	CLFfixval=1.0 ! default value
441	call getin("CLFfixval",CLFfixval)
442	write(,) "CLFfixval=",CLFfixval
443	! sub-grid cloud fraction: varying
444	write(,) "Use partial nebulosity?"
445	CLFvarying=.false. ! default value
446	call getin("CLFvarying",CLFvarying)
447	write(,)"CLFvarying=",CLFvarying
448
449	!CO2 clouds scheme?
450	write(,) "Compute CO2 clouds (implies microphysical scheme)?"
451	co2clouds=.false. ! default value
452	call getin("co2clouds",co2clouds)
453	write(,) " co2clouds = ",co2clouds
454	!Can water ice particles serve as CCN for CO2clouds
455	write(,) "Use water ice as CO2 clouds CCN ?"
456	co2useh2o=.false. ! default value
457	call getin("co2useh2o",co2useh2o)
458	write(,) " co2useh2o = ",co2useh2o
459	!Do we allow a supply of meteoritic paricles to serve as CO2 ice CCN?
460	write(,) "Supply meteoritic particle for CO2 clouds ?"
461	meteo_flux=.false. !Default value
462	call getin("meteo_flux",meteo_flux)
463	write(,) " meteo_flux = ",meteo_flux
464	!Do we allow a sub-grid temperature distribution for the CO2 microphysics
465	write(,) "sub-grid temperature distribution for CO2 clouds?"
466	CLFvaryingCO2=.false. !Default value
467	call getin("CLFvaryingCO2",CLFvaryingCO2)
468	write(,) " CLFvaryingCO2 = ",CLFvaryingCO2
469	!Amplitude of the sub-grid temperature distribution for the CO2 microphysics
470	write(,) "sub-grid temperature amplitude for CO2 clouds?"
471	spantCO2=0 !Default value
472	call getin("spantCO2",spantCO2)
473	write(,) " spantCO2 = ",spantCO2
474
475	! thermal inertia feedback
476	write(,) "Activate the thermal inertia feedback ?"
477	tifeedback=.false. ! default value
478	call getin("tifeedback",tifeedback)
479	write(,) " tifeedback = ",tifeedback
480
481	! Test of incompatibility:
482
483	if (tifeedback.and..not.water) then
484	print*,'if tifeedback is used,'
485	print*,'water should be used too'
486	stop
487	endif
488
489	if (tifeedback.and..not.callsoil) then
490	print*,'if tifeedback is used,'
491	print*,'callsoil should be used too'
492	stop
493	endif
494
495	if (activice.and..not.water) then
496	print*,'if activice is used, water should be used too'
497	stop
498	endif
499
500	if (water.and..not.tracer) then
501	print*,'if water is used, tracer should be used too'
502	stop
503	endif
504
505	! water ice clouds effective variance distribution for sedimentaion
506	write(,) "Sed effective variance for water ice clouds ?"
507	nuice_sed=0.45
508	call getin("nuice_sed",nuice_sed)
509	write(,) "water_param nueff Sedimentation:", nuice_sed
510
511	write(,) "Sed effective variance for CO2 clouds ?"
512	nuiceco2_sed=0.45
513	call getin("nuiceco2_sed",nuiceco2_sed)
514	write(,) "CO2 nueff Sedimentation:", nuiceco2_sed
515
516	write(,) "REF effective variance for CO2 clouds ?"
517	nuiceco2_ref=0.45
518	call getin("nuiceco2_ref",nuiceco2_ref)
519	write(,) "CO2 nueff Sedimentation:", nuiceco2_ref
520
521	write(,) "REF effective variance for water clouds ?"
522	nuice_ref=0.45
523	call getin("nuice_ref",nuice_ref)
524	write(,) "CO2 nueff Sedimentation:", nuice_ref
525
526
527	! ccn factor if no scavenging
528	write(,) "water param CCN reduc. factor ?"
529	ccn_factor = 4.5
530	call getin("ccn_factor",ccn_factor)
531	write(,)" ccn_factor = ",ccn_factor
532	write(,)"Careful: only used when microphys=F, otherwise"
533	write(,)"the contact parameter is used instead;"
534
535	! microphys
536	write(,)"Microphysical scheme for water-ice clouds?"
537	microphys=.false. ! default value
538	call getin("microphys",microphys)
539	write(,)" microphys = ",microphys
540
541	! supersat
542	write(,)"Allow super-saturation of water vapor?"
543	supersat=.true. ! default value
544	call getin("supersat",supersat)
545	write(,)"supersat = ",supersat
546
547	! microphysical parameter contact
548	write(,) "water contact parameter ?"
549	mteta = 0.95
550	call getin("mteta",mteta)
551	write(,) "mteta = ", mteta
552
553	! scavenging
554	write(,)"Dust scavenging by H2O/CO2 snowfall ?"
555	scavenging=.false. ! default value
556	call getin("scavenging",scavenging)
557	write(,)" scavenging = ",scavenging
558
559
560	! Test of incompatibility:
561	! if scavenging is used, then dustbin should be > 0
562
563	if ((microphys.and..not.doubleq).or.
564	& (microphys.and..not.water)) then
565	print*,'if microphys is used, then doubleq,'
566	print*,'and water must be used!'
567	stop
568	endif
569	if (microphys.and..not.scavenging) then
570	print*,''
571	print*,'----------------WARNING-----------------'
572	print*,'microphys is used without scavenging !!!'
573	print*,'----------------WARNING-----------------'
574	print*,''
575	endif
576
577	if ((scavenging.and..not.microphys).or.
578	& (scavenging.and.(dustbin.lt.1)))then
579	print*,'if scavenging is used, then microphys'
580	print*,'must be used!'
581	stop
582	endif
583
584	! Test of incompatibility:
585
586	write(,) "Permanent water caps at poles ?",
587	& " .true. is RECOMMENDED"
588	write(,) "(with .true., North cap is a source of water ",
589	& "and South pole is a cold trap)"
590	caps=.true. ! default value
591	call getin("caps",caps)
592	write(,) " caps = ",caps
593
594	! albedo_h2o_ice
595	write(,) "water ice albedo ?"
596	albedo_h2o_ice=0.45
597	call getin("albedo_h2o_ice",albedo_h2o_ice)
598	write(,) " albedo_h2o_ice = ",albedo_h2o_ice
599	! inert_h2o_ice
600	write(,) "water ice thermal inertia ?"
601	inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
602	call getin("inert_h2o_ice",inert_h2o_ice)
603	write(,) " inert_h2o_ice = ",inert_h2o_ice
604	! frost_albedo_threshold
605	write(,) "frost thickness threshold for albedo ?"
606	frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2)
607	call getin("frost_albedo_threshold",
608	& frost_albedo_threshold)
609	write(,) " frost_albedo_threshold = ",
610	& frost_albedo_threshold
611
612	! call Titus crocus line -- DEFAULT IS NONE
613	write(,) "Titus crocus line ?"
614	tituscap=.false. ! default value
615	call getin("tituscap",tituscap)
616	write(,) "tituscap",tituscap
617
618
619	write(,) "photochemistry: include chemical species"
620	photochem=.false. ! default value
621	call getin("photochem",photochem)
622	write(,) " photochem = ",photochem
623
624
625	! SCATTERERS
626	write(,) "how many scatterers?"
627	naerkind=1 ! default value
628	call getin("naerkind",naerkind)
629	write(,)" naerkind = ",naerkind
630
631	! Test of incompatibility
632	c Logical tests for radiatively active water-ice clouds:
633	IF ( (activice.AND.(.NOT.water)).OR.
634	& (activice.AND.(naerkind.LT.2)) ) THEN
635	WRITE(,) 'If activice is TRUE, water has to be set'
636	WRITE(,) 'to TRUE, and "naerkind" must be at least'
637	WRITE(,) 'equal to 2.'
638	CALL ABORT
639	ELSE IF ( (.NOT.activice).AND.(naerkind.GT.1) ) THEN
640	WRITE(,) 'naerkind is greater than unity, but'
641	WRITE(,) 'activice has not been set to .true.'
642	WRITE(,) 'in callphys.def; this is not logical!'
643	CALL ABORT
644	ENDIF
645
646	!------------------------------------------
647	!------------------------------------------
648	! once naerkind is known allocate arrays
649	! -- we do it here and not in phys_var_init
650	! -- because we need to know naerkind first
651	CALL ini_scatterers(ngrid,nlayer)
652	!------------------------------------------
653	!------------------------------------------
654
655
656	c Please name the different scatterers here ----------------
657	name_iaer(1) = "dust_conrath" !! default choice is good old Conrath profile
658	IF (doubleq.AND.active) name_iaer(1) = "dust_doubleq" !! two-moment scheme
659	if (nq.gt.1) then
660	! trick to avoid problems compiling with 1 tracer
661	! and picky compilers who know name_iaer(2) is out of bounds
662	j=2
663	IF (water.AND.activice) name_iaer(j) = "h2o_ice" !! radiatively-active clouds
664	IF (submicron.AND.active) name_iaer(j) = "dust_submicron" !! JBM experimental stuff
665	endif ! of if (nq.gt.1)
666	c ----------------------------------------------------------
667
668	! THERMOSPHERE
669
670	write(,) "call thermosphere ?"
671	callthermos=.false. ! default value
672	call getin("callthermos",callthermos)
673	write(,) " callthermos = ",callthermos
674
675
676	write(,) " water included without cycle ",
677	& "(only if water=.false.)"
678	thermoswater=.false. ! default value
679	call getin("thermoswater",thermoswater)
680	write(,) " thermoswater = ",thermoswater
681
682	write(,) "call thermal conduction ?",
683	& " (only if callthermos=.true.)"
684	callconduct=.false. ! default value
685	call getin("callconduct",callconduct)
686	write(,) " callconduct = ",callconduct
687
688	write(,) "call EUV heating ?",
689	& " (only if callthermos=.true.)"
690	calleuv=.false. ! default value
691	call getin("calleuv",calleuv)
692	write(,) " calleuv = ",calleuv
693
694	write(,) "call molecular viscosity ?",
695	& " (only if callthermos=.true.)"
696	callmolvis=.false. ! default value
697	call getin("callmolvis",callmolvis)
698	write(,) " callmolvis = ",callmolvis
699
700	write(,) "call molecular diffusion ?",
701	& " (only if callthermos=.true.)"
702	callmoldiff=.false. ! default value
703	call getin("callmoldiff",callmoldiff)
704	write(,) " callmoldiff = ",callmoldiff
705
706
707	write(,) "call thermospheric photochemistry ?",
708	& " (only if callthermos=.true.)"
709	thermochem=.false. ! default value
710	call getin("thermochem",thermochem)
711	write(,) " thermochem = ",thermochem
712
713	write(,) "Method to include solar variability"
714	write(,) "0-> fixed value of E10.7 (fixed_euv_value); ",
715	& "1-> daily evolution of E10.7 (for given solvaryear)"
716	solvarmod=1
717	call getin("solvarmod",solvarmod)
718	write(,) " solvarmod = ",solvarmod
719
720	write(,) "Fixed euv (for solvarmod==0) 10.7 value?"
721	write(,) " (min=80 , ave=140, max=320)"
722	fixed_euv_value=140 ! default value
723	call getin("fixed_euv_value",fixed_euv_value)
724	write(,) " fixed_euv_value = ",fixed_euv_value
725
726	write(,) "Solar variability as observed for MY: "
727	write(,) "Only if solvarmod=1"
728	solvaryear=24
729	call getin("solvaryear",solvaryear)
730	write(,) " solvaryear = ",solvaryear
731
732	write(,) "UV heating efficiency:",
733	& "measured values between 0.19 and 0.23 (Fox et al. 1996)",
734	& "lower values may be used to compensate low 15 um cooling"
735	euveff=0.21 !default value
736	call getin("euveff",euveff)
737	write(,) " euveff = ", euveff
738
739
740	if (.not.callthermos) then
741	if (thermoswater) then
742	print*,'if thermoswater is set, callthermos must be true'
743	stop
744	endif
745	if (callconduct) then
746	print*,'if callconduct is set, callthermos must be true'
747	stop
748	endif
749	if (calleuv) then
750	print*,'if calleuv is set, callthermos must be true'
751	stop
752	endif
753	if (callmolvis) then
754	print*,'if callmolvis is set, callthermos must be true'
755	stop
756	endif
757	if (callmoldiff) then
758	print*,'if callmoldiff is set, callthermos must be true'
759	stop
760	endif
761	if (thermochem) then
762	print*,'if thermochem is set, callthermos must be true'
763	stop
764	endif
765	endif
766
767	! Test of incompatibility:
768	! if photochem is used, then water should be used too
769
770	if (photochem.and..not.water) then
771	print*,'if photochem is used, water should be used too'
772	stop
773	endif
774
775	! if callthermos is used, then thermoswater should be used too
776	! (if water not used already)
777
778	if (callthermos .and. .not.water) then
779	if (callthermos .and. .not.thermoswater) then
780	print*,'if callthermos is used, water or thermoswater
781	& should be used too'
782	stop
783	endif
784	endif
785
786	PRINT*,'--------------------------------------------'
787	PRINT*
788	PRINT*
789	ELSE
790	write(,)
791	write(,) 'Cannot read file callphys.def. Is it here ?'
792	stop
793	ENDIF
794
795	8000 FORMAT(t5,a12,l8)
796	8001 FORMAT(t5,a12,i8)
797
798	PRINT*
799	PRINT*,'conf_phys: daysec',daysec
800	PRINT*
801	PRINT*,'conf_phys: The radiative transfer is computed:'
802	PRINT*,' each ',iradia,' physical time-step'
803	PRINT,' or each ',iradiadtphys,' seconds'
804	PRINT*
805	! --------------------------------------------------------------
806	! Managing the Longwave radiative transfer
807	! --------------------------------------------------------------
808
809	! In most cases, the run just use the following values :
810	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
811	callemis=.true.
812	! ilwd=10*int(daysec/dtphys) ! bug before 22/10/01
813	ilwd=1
814	ilwn=1 !2
815	ilwb=1 !2
816	linear=.true.
817	ncouche=3
818	alphan=0.4
819	semi=0
820
821	! BUT people working hard on the LW may want to read them in 'radia.def'
822	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
823
824	OPEN(99,file='radia.def',status='old',form='formatted'
825	. ,iostat=ierr)
826	IF(ierr.EQ.0) THEN
827	write(,) 'conf_phys: Reading radia.def !!!'
828	READ(99,fmt='(a)') ch1
829	READ(99,*) callemis
830	WRITE(*,8000) ch1,callemis
831
832	READ(99,fmt='(a)') ch1
833	READ(99,*) iradia
834	WRITE(*,8001) ch1,iradia
835
836	READ(99,fmt='(a)') ch1
837	READ(99,*) ilwd
838	WRITE(*,8001) ch1,ilwd
839
840	READ(99,fmt='(a)') ch1
841	READ(99,*) ilwn
842	WRITE(*,8001) ch1,ilwn
843
844	READ(99,fmt='(a)') ch1
845	READ(99,*) linear
846	WRITE(*,8000) ch1,linear
847
848	READ(99,fmt='(a)') ch1
849	READ(99,*) ncouche
850	WRITE(*,8001) ch1,ncouche
851
852	READ(99,fmt='(a)') ch1
853	READ(99,*) alphan
854	WRITE(,) ch1,alphan
855
856	READ(99,fmt='(a)') ch1
857	READ(99,*) ilwb
858	WRITE(*,8001) ch1,ilwb
859
860
861	READ(99,fmt='(a)') ch1
862	READ(99,'(l1)') callg2d
863	WRITE(*,8000) ch1,callg2d
864
865	READ(99,fmt='(a)') ch1
866	READ(99,*) semi
867	WRITE(,) ch1,semi
868	end if
869	CLOSE(99)
870
871	END

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

Download in other formats:

Original Format