Context Navigation

conf_phys.F @ 1711

Last change on this file since 1711 was 1711, checked in by mvals, 8 years ago

Mars GCM

implementation of a sub-grid water cloud fraction scheme (by A. Pottier): If CLFvarying is set to true in callphys.def, then the sub-grid cloud fraction routine is applied in watercloud.F and aeropacity.F.
accordingly modified files are: watercloud.F, aeropacity.F, callradite.F, conf_phys.F, phyetat0.F, phyredem.F90, physiq_mod.F, testphys1d.F, callkeys.h, newtsart.F, start2archive.F, lect_start_archive.F
added file: tcondwater.F90, used by watercloud.F to calculate the condensation temperature of water
watercloud.F, aeropacity.F, callradite.F are converted to module files as watercloud_mod.F, aeropacity_mod.F, callradite_mod.F

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

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

Download in other formats:

Original Format