Context Navigation

inifis.F @ 706

Last change on this file since 706 was 705, checked in by emillour, 12 years ago

Mars GCM:

Added possibility to run with a varying EUV cycle following real one. The flag solvarmod=1 triggers this behaviour, with companion flag solvaryear=## , where ## is the Mars Year (from 23 to 30). Setting solvarmod=0 reverts to 'old' behaviour, where there is a constant EUV forcing throughout the run, set by the "solarcondate" flag.
Needs corresponding input data files ("param_v6" subdirectory of "EUV" subdirectory in "datadir").
Added files jthermcalc_e107.F and param_read_e107.F in "aeronomars", modified files euvheat.F90, hrtherm.F, chemthermos.F90, param_v4.h and param_read.F in "aeronomars" and files inifis.F, physiq.F and callkeys.h in "phymars".

FGG

File size: 27.1 KB

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

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

Download in other formats:

Original Format