Context Navigation

conf_phys.F @ 1544

Last change on this file since 1544 was 1525, checked in by emillour, 9 years ago
All GCMs: More on enforcing dynamics/physics separation: get rid of references to "control_mod" from physics packages. EM
File size: 28.3 KB

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

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

Download in other formats:

Original Format