Context Navigation

conf_phys.F @ 1913

Last change on this file since 1913 was 1818, checked in by jaudouard, 8 years ago
Final update from J.A. about the CO2 clouds scheme for the LMDZ.MARS GCM
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 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	!Do you want to filter the sub-grid T distribution by a Saturation index?
475	write(,) "filter sub-grid temperature by Saturation index?"
476	satindexco2=.true.
477	call getin("satindexco2",satindexco2)
478	write(,) " satindexco2 = ",satindexco2
479
480
481	! thermal inertia feedback
482	write(,) "Activate the thermal inertia feedback ?"
483	tifeedback=.false. ! default value
484	call getin("tifeedback",tifeedback)
485	write(,) " tifeedback = ",tifeedback
486
487	! Test of incompatibility:
488
489	if (tifeedback.and..not.water) then
490	print*,'if tifeedback is used,'
491	print*,'water should be used too'
492	stop
493	endif
494
495	if (tifeedback.and..not.callsoil) then
496	print*,'if tifeedback is used,'
497	print*,'callsoil should be used too'
498	stop
499	endif
500
501	if (activice.and..not.water) then
502	print*,'if activice is used, water should be used too'
503	stop
504	endif
505
506	if (water.and..not.tracer) then
507	print*,'if water is used, tracer should be used too'
508	stop
509	endif
510
511	! water ice clouds effective variance distribution for sedimentaion
512	write(,) "Sed effective variance for water ice clouds ?"
513	nuice_sed=0.45
514	call getin("nuice_sed",nuice_sed)
515	write(,) "water_param nueff Sedimentation:", nuice_sed
516
517	write(,) "Sed effective variance for CO2 clouds ?"
518	nuiceco2_sed=0.45
519	call getin("nuiceco2_sed",nuiceco2_sed)
520	write(,) "CO2 nueff Sedimentation:", nuiceco2_sed
521
522	write(,) "REF effective variance for CO2 clouds ?"
523	nuiceco2_ref=0.45
524	call getin("nuiceco2_ref",nuiceco2_ref)
525	write(,) "CO2 nueff Sedimentation:", nuiceco2_ref
526
527	write(,) "REF effective variance for water clouds ?"
528	nuice_ref=0.45
529	call getin("nuice_ref",nuice_ref)
530	write(,) "CO2 nueff Sedimentation:", nuice_ref
531
532
533	! ccn factor if no scavenging
534	write(,) "water param CCN reduc. factor ?"
535	ccn_factor = 4.5
536	call getin("ccn_factor",ccn_factor)
537	write(,)" ccn_factor = ",ccn_factor
538	write(,)"Careful: only used when microphys=F, otherwise"
539	write(,)"the contact parameter is used instead;"
540
541	! microphys
542	write(,)"Microphysical scheme for water-ice clouds?"
543	microphys=.false. ! default value
544	call getin("microphys",microphys)
545	write(,)" microphys = ",microphys
546
547	! supersat
548	write(,)"Allow super-saturation of water vapor?"
549	supersat=.true. ! default value
550	call getin("supersat",supersat)
551	write(,)"supersat = ",supersat
552
553	! microphysical parameter contact
554	write(,) "water contact parameter ?"
555	mteta = 0.95
556	call getin("mteta",mteta)
557	write(,) "mteta = ", mteta
558
559	! scavenging
560	write(,)"Dust scavenging by H2O/CO2 snowfall ?"
561	scavenging=.false. ! default value
562	call getin("scavenging",scavenging)
563	write(,)" scavenging = ",scavenging
564
565
566	! Test of incompatibility:
567	! if scavenging is used, then dustbin should be > 0
568
569	if ((microphys.and..not.doubleq).or.
570	& (microphys.and..not.water)) then
571	print*,'if microphys is used, then doubleq,'
572	print*,'and water must be used!'
573	stop
574	endif
575	if (microphys.and..not.scavenging) then
576	print*,''
577	print*,'----------------WARNING-----------------'
578	print*,'microphys is used without scavenging !!!'
579	print*,'----------------WARNING-----------------'
580	print*,''
581	endif
582
583	if ((scavenging.and..not.microphys).or.
584	& (scavenging.and.(dustbin.lt.1)))then
585	print*,'if scavenging is used, then microphys'
586	print*,'must be used!'
587	stop
588	endif
589
590	! Test of incompatibility:
591
592	write(,) "Permanent water caps at poles ?",
593	& " .true. is RECOMMENDED"
594	write(,) "(with .true., North cap is a source of water ",
595	& "and South pole is a cold trap)"
596	caps=.true. ! default value
597	call getin("caps",caps)
598	write(,) " caps = ",caps
599
600	! albedo_h2o_ice
601	write(,) "water ice albedo ?"
602	albedo_h2o_ice=0.45
603	call getin("albedo_h2o_ice",albedo_h2o_ice)
604	write(,) " albedo_h2o_ice = ",albedo_h2o_ice
605	! inert_h2o_ice
606	write(,) "water ice thermal inertia ?"
607	inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
608	call getin("inert_h2o_ice",inert_h2o_ice)
609	write(,) " inert_h2o_ice = ",inert_h2o_ice
610	! frost_albedo_threshold
611	write(,) "frost thickness threshold for albedo ?"
612	frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2)
613	call getin("frost_albedo_threshold",
614	& frost_albedo_threshold)
615	write(,) " frost_albedo_threshold = ",
616	& frost_albedo_threshold
617
618	! call Titus crocus line -- DEFAULT IS NONE
619	write(,) "Titus crocus line ?"
620	tituscap=.false. ! default value
621	call getin("tituscap",tituscap)
622	write(,) "tituscap",tituscap
623
624
625	write(,) "photochemistry: include chemical species"
626	photochem=.false. ! default value
627	call getin("photochem",photochem)
628	write(,) " photochem = ",photochem
629
630
631	! SCATTERERS
632	write(,) "how many scatterers?"
633	naerkind=1 ! default value
634	call getin("naerkind",naerkind)
635	write(,)" naerkind = ",naerkind
636
637	! Test of incompatibility
638	c Logical tests for radiatively active water-ice clouds:
639	IF ( (activice.AND.(.NOT.water)).OR.
640	& (activice.AND.(naerkind.LT.2)) ) THEN
641	WRITE(,) 'If activice is TRUE, water has to be set'
642	WRITE(,) 'to TRUE, and "naerkind" must be at least'
643	WRITE(,) 'equal to 2.'
644	CALL ABORT
645	ELSE IF ( (.NOT.activice).AND.(naerkind.GT.1) ) THEN
646	WRITE(,) 'naerkind is greater than unity, but'
647	WRITE(,) 'activice has not been set to .true.'
648	WRITE(,) 'in callphys.def; this is not logical!'
649	CALL ABORT
650	ENDIF
651
652	!------------------------------------------
653	!------------------------------------------
654	! once naerkind is known allocate arrays
655	! -- we do it here and not in phys_var_init
656	! -- because we need to know naerkind first
657	CALL ini_scatterers(ngrid,nlayer)
658	!------------------------------------------
659	!------------------------------------------
660
661
662	c Please name the different scatterers here ----------------
663	name_iaer(1) = "dust_conrath" !! default choice is good old Conrath profile
664	IF (doubleq.AND.active) name_iaer(1) = "dust_doubleq" !! two-moment scheme
665	if (nq.gt.1) then
666	! trick to avoid problems compiling with 1 tracer
667	! and picky compilers who know name_iaer(2) is out of bounds
668	j=2
669	IF (water.AND.activice) name_iaer(j) = "h2o_ice" !! radiatively-active clouds
670	IF (submicron.AND.active) name_iaer(j) = "dust_submicron" !! JBM experimental stuff
671	endif ! of if (nq.gt.1)
672	c ----------------------------------------------------------
673
674	! THERMOSPHERE
675
676	write(,) "call thermosphere ?"
677	callthermos=.false. ! default value
678	call getin("callthermos",callthermos)
679	write(,) " callthermos = ",callthermos
680
681
682	write(,) " water included without cycle ",
683	& "(only if water=.false.)"
684	thermoswater=.false. ! default value
685	call getin("thermoswater",thermoswater)
686	write(,) " thermoswater = ",thermoswater
687
688	write(,) "call thermal conduction ?",
689	& " (only if callthermos=.true.)"
690	callconduct=.false. ! default value
691	call getin("callconduct",callconduct)
692	write(,) " callconduct = ",callconduct
693
694	write(,) "call EUV heating ?",
695	& " (only if callthermos=.true.)"
696	calleuv=.false. ! default value
697	call getin("calleuv",calleuv)
698	write(,) " calleuv = ",calleuv
699
700	write(,) "call molecular viscosity ?",
701	& " (only if callthermos=.true.)"
702	callmolvis=.false. ! default value
703	call getin("callmolvis",callmolvis)
704	write(,) " callmolvis = ",callmolvis
705
706	write(,) "call molecular diffusion ?",
707	& " (only if callthermos=.true.)"
708	callmoldiff=.false. ! default value
709	call getin("callmoldiff",callmoldiff)
710	write(,) " callmoldiff = ",callmoldiff
711
712
713	write(,) "call thermospheric photochemistry ?",
714	& " (only if callthermos=.true.)"
715	thermochem=.false. ! default value
716	call getin("thermochem",thermochem)
717	write(,) " thermochem = ",thermochem
718
719	write(,) "Method to include solar variability"
720	write(,) "0-> fixed value of E10.7 (fixed_euv_value); ",
721	& "1-> daily evolution of E10.7 (for given solvaryear)"
722	solvarmod=1
723	call getin("solvarmod",solvarmod)
724	write(,) " solvarmod = ",solvarmod
725
726	write(,) "Fixed euv (for solvarmod==0) 10.7 value?"
727	write(,) " (min=80 , ave=140, max=320)"
728	fixed_euv_value=140 ! default value
729	call getin("fixed_euv_value",fixed_euv_value)
730	write(,) " fixed_euv_value = ",fixed_euv_value
731
732	write(,) "Solar variability as observed for MY: "
733	write(,) "Only if solvarmod=1"
734	solvaryear=24
735	call getin("solvaryear",solvaryear)
736	write(,) " solvaryear = ",solvaryear
737
738	write(,) "UV heating efficiency:",
739	& "measured values between 0.19 and 0.23 (Fox et al. 1996)",
740	& "lower values may be used to compensate low 15 um cooling"
741	euveff=0.21 !default value
742	call getin("euveff",euveff)
743	write(,) " euveff = ", euveff
744
745
746	if (.not.callthermos) then
747	if (thermoswater) then
748	print*,'if thermoswater is set, callthermos must be true'
749	stop
750	endif
751	if (callconduct) then
752	print*,'if callconduct is set, callthermos must be true'
753	stop
754	endif
755	if (calleuv) then
756	print*,'if calleuv is set, callthermos must be true'
757	stop
758	endif
759	if (callmolvis) then
760	print*,'if callmolvis is set, callthermos must be true'
761	stop
762	endif
763	if (callmoldiff) then
764	print*,'if callmoldiff is set, callthermos must be true'
765	stop
766	endif
767	if (thermochem) then
768	print*,'if thermochem is set, callthermos must be true'
769	stop
770	endif
771	endif
772
773	! Test of incompatibility:
774	! if photochem is used, then water should be used too
775
776	if (photochem.and..not.water) then
777	print*,'if photochem is used, water should be used too'
778	stop
779	endif
780
781	! if callthermos is used, then thermoswater should be used too
782	! (if water not used already)
783
784	if (callthermos .and. .not.water) then
785	if (callthermos .and. .not.thermoswater) then
786	print*,'if callthermos is used, water or thermoswater
787	& should be used too'
788	stop
789	endif
790	endif
791
792	PRINT*,'--------------------------------------------'
793	PRINT*
794	PRINT*
795	ELSE
796	write(,)
797	write(,) 'Cannot read file callphys.def. Is it here ?'
798	stop
799	ENDIF
800
801	8000 FORMAT(t5,a12,l8)
802	8001 FORMAT(t5,a12,i8)
803
804	PRINT*
805	PRINT*,'conf_phys: daysec',daysec
806	PRINT*
807	PRINT*,'conf_phys: The radiative transfer is computed:'
808	PRINT*,' each ',iradia,' physical time-step'
809	PRINT,' or each ',iradiadtphys,' seconds'
810	PRINT*
811	! --------------------------------------------------------------
812	! Managing the Longwave radiative transfer
813	! --------------------------------------------------------------
814
815	! In most cases, the run just use the following values :
816	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
817	callemis=.true.
818	! ilwd=10*int(daysec/dtphys) ! bug before 22/10/01
819	ilwd=1
820	ilwn=1 !2
821	ilwb=1 !2
822	linear=.true.
823	ncouche=3
824	alphan=0.4
825	semi=0
826
827	! BUT people working hard on the LW may want to read them in 'radia.def'
828	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
829
830	OPEN(99,file='radia.def',status='old',form='formatted'
831	. ,iostat=ierr)
832	IF(ierr.EQ.0) THEN
833	write(,) 'conf_phys: Reading radia.def !!!'
834	READ(99,fmt='(a)') ch1
835	READ(99,*) callemis
836	WRITE(*,8000) ch1,callemis
837
838	READ(99,fmt='(a)') ch1
839	READ(99,*) iradia
840	WRITE(*,8001) ch1,iradia
841
842	READ(99,fmt='(a)') ch1
843	READ(99,*) ilwd
844	WRITE(*,8001) ch1,ilwd
845
846	READ(99,fmt='(a)') ch1
847	READ(99,*) ilwn
848	WRITE(*,8001) ch1,ilwn
849
850	READ(99,fmt='(a)') ch1
851	READ(99,*) linear
852	WRITE(*,8000) ch1,linear
853
854	READ(99,fmt='(a)') ch1
855	READ(99,*) ncouche
856	WRITE(*,8001) ch1,ncouche
857
858	READ(99,fmt='(a)') ch1
859	READ(99,*) alphan
860	WRITE(,) ch1,alphan
861
862	READ(99,fmt='(a)') ch1
863	READ(99,*) ilwb
864	WRITE(*,8001) ch1,ilwb
865
866
867	READ(99,fmt='(a)') ch1
868	READ(99,'(l1)') callg2d
869	WRITE(*,8000) ch1,callg2d
870
871	READ(99,fmt='(a)') ch1
872	READ(99,*) semi
873	WRITE(,) ch1,semi
874	end if
875	CLOSE(99)
876
877	END

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

Download in other formats:

Original Format