Context Navigation

conf_phys.F @ 1617

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

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

Download in other formats:

Original Format