Context Navigation

inifis.F @ 1212

Last change on this file since 1212 was 1212, checked in by aslmd, 11 years ago

LMDZ.MARS + MESOSCALE

A quite major commit, at least for MESOSCALE.
In a word: any ngrid deserves to be free.

no need to recompile when changing number of horizontal grid points or number of processors
latest version of LMDZ.MARS physics can be used
WARNING! Nesting is still yet to be fixed (since r1027)

Also some small bug fixes to LMDZ.MARS.

Changes in LMDZ.MARS

--> fixed a potential bug in thermal plume model because zlmax was computed both in thermcell_main_mars and calltherm_interface... so made it an OUT argument of calltherm_interface. also: changed the name to limz. and added precompiling flags to avoid the use of planetwide in MESOSCALE. in MESOSCALE we just go high enough (nlayer-5) and do not care about computational cost (although we certainly gain from not using MAXVAL).
--> moved allocations upward in inifis. does not change anything for GCM, but make MESOSCALE modifications simpler, and overall make inifis better organized: first allocations, then reading callphys.def file.
--> added precompiling flags around lines that are both useless for MESOSCALE (notably I/O) and recently adapted to parallel computations in the GCM
--> tidied up what is MESOSCALE vs. GCM in surfini

Changes in MESOSCALE

--> changed makemeso to allow dynamically set nx ny nprocs
--> changed makemeso to remove links to Fortran code adapted to parallel GCM and useless for mesoscale
--> changed ngridmx to ngrid in inifis includes

File size: 27.9 KB

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

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

Download in other formats:

Original Format