Context Navigation

source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/phys/inifis.f @ 224

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago

File size: 30.6 KB

Line
1	SUBROUTINE inifis(ngrid,nlayer,nq,
2	$ day_ini,pdaysec,ptimestep,
3	$ plat,plon,parea,
4	$ prad,pg,pr,pcpp)
5
6	use radinc_h, only : naerkind
7	use datafile_mod, only: datadir
8	use comdiurn_h, only: sinlat, coslat, sinlon, coslon
9	use comgeomfi_h, only: long, lati, area, totarea, totarea_planet
10	use comsoil_h, only: ini_comsoil_h
11	use control_mod, only: ecritphy
12	use planetwide_mod, only: planetwide_sumval
13
14	!=======================================================================
15	!
16	! purpose:
17	! -------
18	!
19	! Initialisation for the physical parametrisations of the LMD
20	! Generic Model.
21	!
22	! author: Frederic Hourdin 15 / 10 /93
23	! -------
24	! modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
25	! Ehouarn Millour (oct. 2008) tracers are now identified
26	! by their names and may not be contiguously
27	! stored in the q(:,:,:,:) array
28	! E.M. (june 2009) use getin routine to load parameters
29	!
30	!
31	! arguments:
32	! ----------
33	!
34	! input:
35	! ------
36	!
37	! ngrid Size of the horizontal grid.
38	! All internal loops are performed on that grid.
39	! nlayer Number of vertical layers.
40	! pdayref Day of reference for the simulation
41	! pday Number of days counted from the North. Spring
42	! equinoxe.
43	!
44	!=======================================================================
45	!
46	!-----------------------------------------------------------------------
47	! declarations:
48	! -------------
49	use datafile_mod, only: datadir
50	! to use 'getin'
51	USE ioipsl_getincom
52	IMPLICIT NONE
53	!-----------------------------------------------------------------------
54	! INCLUDE 'dimensions.h'
55	!
56	! dimensions.h contient les dimensions du modele
57	! ndm est tel que iim=2**ndm
58	!-----------------------------------------------------------------------
59
60	INTEGER iim,jjm,llm,ndm
61
62	PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
63
64	!-----------------------------------------------------------------------
65	!-----------------------------------------------------------------------
66	! INCLUDE 'dimphys.h'
67
68	! ngridmx : number of horizontal grid points
69	! note: the -1/jjm term will be 0; unless jj=1
70	integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)
71	! nlayermx : number of atmospheric layers
72	integer, parameter :: nlayermx = llm
73	! nsoilmx : number of subterranean layers ! nsoilmx is now in comsoil_h
74	!integer, parameter :: nsoilmx = 4 ! for a test
75	!integer, parameter :: nsoilmx = 18 ! for z1=0.0002 m, depth = 18 m => mars case
76	!integer, parameter :: nsoilmx = 13 ! for z1=0.03 m, depth = 104.8 m => earth case
77	!-----------------------------------------------------------------------
78
79	!-----------------------------------------------------------------------
80	! INCLUDE planet.h
81
82	COMMON/planet/apoastr,periastr,year_day,peri_day, &
83	& obliquit,nres, &
84	& z0,lmixmin,emin_turb,coefvis,coefir, &
85	& timeperi,e_elips,p_elips
86
87	real apoastr,periastr,year_day,peri_day, &
88	& obliquit,nres, &
89	& z0,lmixmin,emin_turb,coefvis,coefir, &
90	& timeperi,e_elips,p_elips
91
92
93	!-----------------------------------------------------------------------
94	!-----------------------------------------------------------------------
95	! INCLUDE "comcstfi.h"
96
97	common/comcstfi/pi,rad,g,r,cpp,rcp,dtphys,daysec,mugaz,omeg
98	common/comcstfi/avocado!,molrad,visc
99
100	real pi,rad,g,r,cpp,rcp,dtphys,daysec,mugaz,omeg
101	real avocado!,molrad,visc
102
103	!
104	! For Fortran 77/Fortran 90 compliance always use line continuation
105	! symbols '&' in columns 73 and 6
106	!
107	! Group commons according to their type for minimal performance impact
108
109	COMMON/callkeys_l/callrad,corrk,calldifv,UseTurbDiff,calladj &
110	& , co2cond,callsoil &
111	& , season,diurnal,tlocked,rings_shadow,lwrite &
112	& , callstats,calleofdump &
113	& , enertest &
114	& , callgasvis,continuum,H2Ocont_simple,graybody &
115	& , radfixed &
116	& , meanOLR, specOLR &
117	& , kastprof &
118	& , nosurf, oblate &
119	& , newtonian, testradtimes &
120	& , check_cpp_match, force_cpp &
121	& , rayleigh &
122	& , stelbbody &
123	& , nearco2cond &
124	& , tracer, mass_redistrib, varactive, varfixed &
125	& , sedimentation,water,watercond,waterrain &
126	& , aeroco2,aeroh2o,aeroh2so4,aeroback2lay &
127	& , aerofixco2,aerofixh2o &
128	& , hydrology, sourceevol &
129	& , CLFvarying &
130	& , strictboundcorrk &
131	& , ok_slab_ocean &
132	& , ok_slab_sic &
133	& , ok_slab_heat_transp
134
135
136	COMMON/callkeys_i/iaervar,iddist,iradia,startype
137
138	COMMON/callkeys_r/topdustref,Nmix_co2,dusttau,Fat1AU,stelTbb, &
139	& Tstrat,tplanet,obs_tau_col_tropo, &
140	& obs_tau_col_strato,pres_bottom_tropo, &
141	& pres_top_tropo,pres_bottom_strato, &
142	& pres_top_strato,size_tropo,size_strato,satval, &
143	& CLFfixval,n2mixratio,co2supsat,pceil,albedosnow,&
144	& maxicethick,Tsaldiff,tau_relax,cloudlvl, &
145	& icetstep,intheat,flatten,Rmean,J2,MassPlanet
146
147	logical callrad,corrk,calldifv,UseTurbDiff &
148	& , calladj,co2cond,callsoil &
149	& , season,diurnal,tlocked,rings_shadow,lwrite &
150	& , callstats,calleofdump &
151	& , callgasvis,continuum,H2Ocont_simple,graybody &
152	& , strictboundcorrk
153
154	logical enertest
155	logical nonideal
156	logical meanOLR
157	logical specOLR
158	logical kastprof
159	logical newtonian
160	logical check_cpp_match
161	logical force_cpp
162	logical testradtimes
163	logical rayleigh
164	logical stelbbody
165	logical ozone
166	logical nearco2cond
167	logical tracer
168	logical mass_redistrib
169	logical varactive
170	logical varfixed
171	logical radfixed
172	logical sedimentation
173	logical water,watercond,waterrain
174	logical aeroco2,aeroh2o,aeroh2so4,aeroback2lay
175	logical aerofixco2,aerofixh2o
176	logical hydrology
177	logical sourceevol
178	logical CLFvarying
179	logical nosurf
180	logical oblate
181	logical ok_slab_ocean
182	logical ok_slab_sic
183	logical ok_slab_heat_transp
184
185	integer iddist
186	integer iaervar
187	integer iradia
188	integer startype
189
190	real topdustref
191	real Nmix_co2
192	real dusttau
193	real Fat1AU
194	real stelTbb
195	real Tstrat
196	real tplanet
197	real obs_tau_col_tropo
198	real obs_tau_col_strato
199	real pres_bottom_tropo
200	real pres_top_tropo
201	real pres_bottom_strato
202	real pres_top_strato
203	real size_tropo
204	real size_strato
205	real satval
206	real CLFfixval
207	real n2mixratio
208	real co2supsat
209	real pceil
210	real albedosnow
211	real maxicethick
212	real Tsaldiff
213	real tau_relax
214	real cloudlvl
215	real icetstep
216	real intheat
217	real flatten
218	real Rmean
219	real J2
220	real MassPlanet
221
222
223
224	REAL,INTENT(IN) :: prad,pg,pr,pcpp,pdaysec,ptimestep
225
226	INTEGER,INTENT(IN) :: ngrid,nlayer,nq
227	REAL,INTENT(IN) :: plat(ngrid),plon(ngrid),parea(ngrid)
228	integer day_ini
229	INTEGER ig,ierr
230
231	EXTERNAL iniorbit,orbite
232	EXTERNAL SSUM
233	REAL SSUM
234
235	CHARACTER ch1*12
236	CHARACTER ch80*80
237
238	logical chem, h2o
239	logical :: parameter, doubleq=.false.
240
241	real psurf,pN2 ! added by RW for Gliese 581d N2+CO2
242
243	rad=prad
244	daysec=pdaysec
245	dtphys=ptimestep
246	cpp=pcpp
247	g=pg
248	r=pr
249	rcp=r/cpp
250
251	avocado = 6.02214179e23 ! added by RW
252
253	! --------------------------------------------------------
254	! The usual Tests
255	! --------------
256	IF (nlayer.NE.nlayermx) THEN
257	PRINT*,'STOP in inifis'
258	PRINT*,'Probleme de dimensions :'
259	PRINT*,'nlayer = ',nlayer
260	PRINT*,'nlayermx = ',nlayermx
261	STOP
262	ENDIF
263
264	! read in 'ecritphy' (frequency of calls to physics, in dynamical steps)
265	! (also done in dyn3d/defrun_new but not in LMDZ.COMMON)
266	call getin("ecritphy",ecritphy)
267
268	! --------------------------------------------------------------
269	! Reading the "callphys.def" file controlling some key options
270	! --------------------------------------------------------------
271
272	! check that 'callphys.def' file is around
273	OPEN(99,file='callphys.def',status='old',form='formatted'
274	& ,iostat=ierr)
275	CLOSE(99)
276
277	IF(ierr.EQ.0) THEN
278	PRINT*
279	PRINT*
280	PRINT*,'--------------------------------------------'
281	PRINT*,' inifis: Parametres pour la physique (callphys.def)'
282	PRINT*,'--------------------------------------------'
283
284	write(,) "Directory where external input files are:"
285	! default 'datadir' is set in "datadir_mod"
286	call getin("datadir",datadir) ! default path
287	write(,) " datadir = ",trim(datadir)
288
289	write(,) "Run with or without tracer transport ?"
290	tracer=.false. ! default value
291	call getin("tracer",tracer)
292	write(,) " tracer = ",tracer
293
294	write(,) "Run with or without atm mass update ",
295	& " due to tracer evaporation/condensation?"
296	mass_redistrib=.false. ! default value
297	call getin("mass_redistrib",mass_redistrib)
298	write(,) " mass_redistrib = ",mass_redistrib
299
300	write(,) "Diurnal cycle ?"
301	write(,) "(if diurnal=false, diurnal averaged solar heating)"
302	diurnal=.true. ! default value
303	call getin("diurnal",diurnal)
304	write(,) " diurnal = ",diurnal
305
306	write(,) "Seasonal cycle ?"
307	write(,) "(if season=false, Ls stays constant, to value ",
308	& "set in 'start'"
309	season=.true. ! default value
310	call getin("season",season)
311	write(,) " season = ",season
312
313	write(,) "Tidally resonant rotation ?"
314	tlocked=.false. ! default value
315	call getin("tlocked",tlocked)
316	write(,) "tlocked = ",tlocked
317
318	write(,) "Saturn ring shadowing ?"
319	rings_shadow = .false.
320	call getin("rings_shadow", rings_shadow)
321	write(,) "rings_shadow = ", rings_shadow
322
323	write(,) "Compute latitude-dependent gravity field?"
324	oblate = .false.
325	call getin("oblate", oblate)
326	write(,) "oblate = ", oblate
327
328	write(,) "Flattening of the planet (a-b)/a "
329	flatten = 0.0
330	call getin("flatten", flatten)
331	write(,) "flatten = ", flatten
332
333
334	write(,) "Needed if oblate=.true.: J2"
335	J2 = 0.0
336	call getin("J2", J2)
337	write(,) "J2 = ", J2
338
339	write(,) "Needed if oblate=.true.: Planet mass (*1e24 kg)"
340	MassPlanet = 0.0
341	call getin("MassPlanet", MassPlanet)
342	write(,) "MassPlanet = ", MassPlanet
343
344	write(,) "Needed if oblate=.true.: Planet mean radius (m)"
345	Rmean = 0.0
346	call getin("Rmean", Rmean)
347	write(,) "Rmean = ", Rmean
348
349	! Test of incompatibility:
350	! if tlocked, then diurnal should be false
351	if (tlocked.and.diurnal) then
352	print*,'If diurnal=true, we should turn off tlocked.'
353	stop
354	endif
355
356	write(,) "Tidal resonance ratio ?"
357	nres=0 ! default value
358	call getin("nres",nres)
359	write(,) "nres = ",nres
360
361	write(,) "Write some extra output to the screen ?"
362	lwrite=.false. ! default value
363	call getin("lwrite",lwrite)
364	write(,) " lwrite = ",lwrite
365
366	write(,) "Save statistics in file stats.nc ?"
367	callstats=.true. ! default value
368	call getin("callstats",callstats)
369	write(,) " callstats = ",callstats
370
371	write(,) "Test energy conservation of model physics ?"
372	enertest=.false. ! default value
373	call getin("enertest",enertest)
374	write(,) " enertest = ",enertest
375
376	write(,) "Check to see if cpp values used match gases.def ?"
377	check_cpp_match=.true. ! default value
378	call getin("check_cpp_match",check_cpp_match)
379	write(,) " check_cpp_match = ",check_cpp_match
380
381	write(,) "call radiative transfer ?"
382	callrad=.true. ! default value
383	call getin("callrad",callrad)
384	write(,) " callrad = ",callrad
385
386	write(,) "call correlated-k radiative transfer ?"
387	corrk=.true. ! default value
388	call getin("corrk",corrk)
389	write(,) " corrk = ",corrk
390
391	write(,) "prohibit calculations outside corrk T grid?"
392	strictboundcorrk=.true. ! default value
393	call getin("strictboundcorrk",strictboundcorrk)
394	write(,) "strictboundcorrk = ",strictboundcorrk
395
396	write(,) "call gaseous absorption in the visible bands?",
397	& "(matters only if callrad=T)"
398	callgasvis=.false. ! default value
399	call getin("callgasvis",callgasvis)
400	write(,) " callgasvis = ",callgasvis
401
402	write(,) "call continuum opacities in radiative transfer ?",
403	& "(matters only if callrad=T)"
404	continuum=.true. ! default value
405	call getin("continuum",continuum)
406	write(,) " continuum = ",continuum
407
408	write(,) "use analytic function for H2O continuum ?"
409	H2Ocont_simple=.false. ! default value
410	call getin("H2Ocont_simple",H2Ocont_simple)
411	write(,) " H2Ocont_simple = ",H2Ocont_simple
412
413	write(,) "call turbulent vertical diffusion ?"
414	calldifv=.true. ! default value
415	call getin("calldifv",calldifv)
416	write(,) " calldifv = ",calldifv
417
418	write(,) "use turbdiff instead of vdifc ?"
419	UseTurbDiff=.true. ! default value
420	call getin("UseTurbDiff",UseTurbDiff)
421	write(,) " UseTurbDiff = ",UseTurbDiff
422
423	write(,) "call convective adjustment ?"
424	calladj=.true. ! default value
425	call getin("calladj",calladj)
426	write(,) " calladj = ",calladj
427
428	write(,) "call CO2 condensation ?"
429	co2cond=.false. ! default value
430	call getin("co2cond",co2cond)
431	write(,) " co2cond = ",co2cond
432	! Test of incompatibility
433	if (co2cond.and.(.not.tracer)) then
434	print*,'We need a CO2 ice tracer to condense CO2'
435	call abort
436	endif
437
438	write(,) "CO2 supersaturation level ?"
439	co2supsat=1.0 ! default value
440	call getin("co2supsat",co2supsat)
441	write(,) " co2supsat = ",co2supsat
442
443	write(,) "Radiative timescale for Newtonian cooling ?"
444	tau_relax=30. ! default value
445	call getin("tau_relax",tau_relax)
446	write(,) " tau_relax = ",tau_relax
447	tau_relax=tau_relax243600 ! convert Earth days --> seconds
448
449	write(,)"call thermal conduction in the soil ?"
450	callsoil=.true. ! default value
451	call getin("callsoil",callsoil)
452	write(,) " callsoil = ",callsoil
453
454	write(,)"Rad transfer is computed every iradia",
455	& " physical timestep"
456	iradia=1 ! default value
457	call getin("iradia",iradia)
458	write(,)" iradia = ",iradia
459
460	write(,)"Rayleigh scattering ?"
461	rayleigh=.false.
462	call getin("rayleigh",rayleigh)
463	write(,)" rayleigh = ",rayleigh
464
465	write(,) "Use blackbody for stellar spectrum ?"
466	stelbbody=.false. ! default value
467	call getin("stelbbody",stelbbody)
468	write(,) " stelbbody = ",stelbbody
469
470	write(,) "Stellar blackbody temperature ?"
471	stelTbb=5800.0 ! default value
472	call getin("stelTbb",stelTbb)
473	write(,) " stelTbb = ",stelTbb
474
475	write(,)"Output mean OLR in 1D?"
476	meanOLR=.false.
477	call getin("meanOLR",meanOLR)
478	write(,)" meanOLR = ",meanOLR
479
480	write(,)"Output spectral OLR in 3D?"
481	specOLR=.false.
482	call getin("specOLR",specOLR)
483	write(,)" specOLR = ",specOLR
484
485	write(,)"Operate in kastprof mode?"
486	kastprof=.false.
487	call getin("kastprof",kastprof)
488	write(,)" kastprof = ",kastprof
489
490	write(,)"Uniform absorption in radiative transfer?"
491	graybody=.false.
492	call getin("graybody",graybody)
493	write(,)" graybody = ",graybody
494
495	! Slab Ocean
496	write(,) "Use slab-ocean ?"
497	ok_slab_ocean=.false. ! default value
498	call getin("ok_slab_ocean",ok_slab_ocean)
499	write(,) "ok_slab_ocean = ",ok_slab_ocean
500
501	write(,) "Use slab-sea-ice ?"
502	ok_slab_sic=.true. ! default value
503	call getin("ok_slab_sic",ok_slab_sic)
504	write(,) "ok_slab_sic = ",ok_slab_sic
505
506	write(,) "Use heat transport for the ocean ?"
507	ok_slab_heat_transp=.true. ! default value
508	call getin("ok_slab_heat_transp",ok_slab_heat_transp)
509	write(,) "ok_slab_heat_transp = ",ok_slab_heat_transp
510
511
512
513	! Test of incompatibility:
514	! if kastprof used, we must be in 1D
515	if (kastprof.and.(ngrid.gt.1)) then
516	print*,'kastprof can only be used in 1D!'
517	call abort
518	endif
519
520	write(,)"Stratospheric temperature for kastprof mode?"
521	Tstrat=167.0
522	call getin("Tstrat",Tstrat)
523	write(,)" Tstrat = ",Tstrat
524
525	write(,)"Remove lower boundary?"
526	nosurf=.false.
527	call getin("nosurf",nosurf)
528	write(,)" nosurf = ",nosurf
529
530	! Tests of incompatibility:
531	if (nosurf.and.callsoil) then
532	print*,'nosurf not compatible with soil scheme!'
533	print*,'... got to make a choice!'
534	call abort
535	endif
536
537	write(,)"Add an internal heat flux?",
538	. "... matters only if callsoil=F"
539	intheat=0.
540	call getin("intheat",intheat)
541	write(,)" intheat = ",intheat
542
543	write(,)"Use Newtonian cooling for radiative transfer?"
544	newtonian=.false.
545	call getin("newtonian",newtonian)
546	write(,)" newtonian = ",newtonian
547
548	! Tests of incompatibility:
549	if (newtonian.and.corrk) then
550	print*,'newtonian not compatible with correlated-k!'
551	call abort
552	endif
553	if (newtonian.and.calladj) then
554	print*,'newtonian not compatible with adjustment!'
555	call abort
556	endif
557	if (newtonian.and.calldifv) then
558	print*,'newtonian not compatible with a boundary layer!'
559	call abort
560	endif
561
562	write(,)"Test physics timescale in 1D?"
563	testradtimes=.false.
564	call getin("testradtimes",testradtimes)
565	write(,)" testradtimes = ",testradtimes
566
567	! Test of incompatibility:
568	! if testradtimes used, we must be in 1D
569	if (testradtimes.and.(ngrid.gt.1)) then
570	print*,'testradtimes can only be used in 1D!'
571	call abort
572	endif
573
574	write(,)"Default planetary temperature?"
575	tplanet=215.0
576	call getin("tplanet",tplanet)
577	write(,)" tplanet = ",tplanet
578
579	write(,)"Which star?"
580	startype=1 ! default value = Sol
581	call getin("startype",startype)
582	write(,)" startype = ",startype
583
584	write(,)"Value of stellar flux at 1 AU?"
585	Fat1AU=1356.0 ! default value = Sol today
586	call getin("Fat1AU",Fat1AU)
587	write(,)" Fat1AU = ",Fat1AU
588
589
590	! TRACERS:
591
592	write(,)"Varying H2O cloud fraction?"
593	CLFvarying=.false. ! default value
594	call getin("CLFvarying",CLFvarying)
595	write(,)" CLFvarying = ",CLFvarying
596
597	write(,)"Value of fixed H2O cloud fraction?"
598	CLFfixval=1.0 ! default value
599	call getin("CLFfixval",CLFfixval)
600	write(,)" CLFfixval = ",CLFfixval
601
602	write(,)"fixed radii for Cloud particles?"
603	radfixed=.false. ! default value
604	call getin("radfixed",radfixed)
605	write(,)" radfixed = ",radfixed
606
607	if(kastprof)then
608	radfixed=.true.
609	endif
610
611	write(,)"Number mixing ratio of CO2 ice particles:"
612	Nmix_co2=1.e6 ! default value
613	call getin("Nmix_co2",Nmix_co2)
614	write(,)" Nmix_co2 = ",Nmix_co2
615
616	! write(,)"Number of radiatively active aerosols:"
617	! naerkind=0. ! default value
618	! call getin("naerkind",naerkind)
619	! write(,)" naerkind = ",naerkind
620
621	write(,)"Opacity of dust (if used):"
622	dusttau=0. ! default value
623	call getin("dusttau",dusttau)
624	write(,)" dusttau = ",dusttau
625
626	write(,)"Radiatively active CO2 aerosols?"
627	aeroco2=.false. ! default value
628	call getin("aeroco2",aeroco2)
629	write(,)" aeroco2 = ",aeroco2
630
631	write(,)"Fixed CO2 aerosol distribution?"
632	aerofixco2=.false. ! default value
633	call getin("aerofixco2",aerofixco2)
634	write(,)" aerofixco2 = ",aerofixco2
635
636	write(,)"Radiatively active water ice?"
637	aeroh2o=.false. ! default value
638	call getin("aeroh2o",aeroh2o)
639	write(,)" aeroh2o = ",aeroh2o
640
641	write(,)"Fixed H2O aerosol distribution?"
642	aerofixh2o=.false. ! default value
643	call getin("aerofixh2o",aerofixh2o)
644	write(,)" aerofixh2o = ",aerofixh2o
645
646	write(,)"Radiatively active sulfuric acid aersols?"
647	aeroh2so4=.false. ! default value
648	call getin("aeroh2so4",aeroh2so4)
649	write(,)" aeroh2so4 = ",aeroh2so4
650
651	!=================================
652
653	write(,)"Radiatively active two-layer aersols?"
654	aeroback2lay=.false. ! default value
655	call getin("aeroback2lay",aeroback2lay)
656	write(,)" aeroback2lay = ",aeroback2lay
657
658	write(,)"TWOLAY AEROSOL: total optical depth ",
659	& "in the tropospheric layer (visible)"
660	obs_tau_col_tropo=8.D0
661	call getin("obs_tau_col_tropo",obs_tau_col_tropo)
662	write(,)" obs_tau_col_tropo = ",obs_tau_col_tropo
663
664	write(,)"TWOLAY AEROSOL: total optical depth ",
665	& "in the stratospheric layer (visible)"
666	obs_tau_col_strato=0.08D0
667	call getin("obs_tau_col_strato",obs_tau_col_strato)
668	write(,)" obs_tau_col_strato = ",obs_tau_col_strato
669
670	write(,)"TWOLAY AEROSOL: pres_bottom_tropo? in pa"
671	pres_bottom_tropo=66000.0
672	call getin("pres_bottom_tropo",pres_bottom_tropo)
673	write(,)" pres_bottom_tropo = ",pres_bottom_tropo
674
675	write(,)"TWOLAY AEROSOL: pres_top_tropo? in pa"
676	pres_top_tropo=18000.0
677	call getin("pres_top_tropo",pres_top_tropo)
678	write(,)" pres_top_tropo = ",pres_top_tropo
679
680	write(,)"TWOLAY AEROSOL: pres_bottom_strato? in pa"
681	pres_bottom_strato=2000.0
682	call getin("pres_bottom_strato",pres_bottom_strato)
683	write(,)" pres_bottom_strato = ",pres_bottom_strato
684
685	write(,)"TWOLAY AEROSOL: pres_top_strato? in pa"
686	pres_top_strato=100.0
687	call getin("pres_top_strato",pres_top_strato)
688	write(,)" pres_top_strato = ",pres_top_strato
689
690	write(,)"TWOLAY AEROSOL: particle size in the ",
691	& "tropospheric layer, in meters"
692	size_tropo=2.e-6
693	call getin("size_tropo",size_tropo)
694	write(,)" size_tropo = ",size_tropo
695
696	write(,)"TWOLAY AEROSOL: particle size in the ",
697	& "stratospheric layer, in meters"
698	size_strato=1.e-7
699	call getin("size_strato",size_strato)
700	write(,)" size_strato = ",size_strato
701
702	!=================================
703
704	write(,)"Cloud pressure level (with kastprof only):"
705	cloudlvl=0. ! default value
706	call getin("cloudlvl",cloudlvl)
707	write(,)" cloudlvl = ",cloudlvl
708
709	write(,)"Is the variable gas species radiatively active?"
710	Tstrat=167.0
711	varactive=.false.
712	call getin("varactive",varactive)
713	write(,)" varactive = ",varactive
714
715	write(,)"Is the variable gas species distribution set?"
716	varfixed=.false.
717	call getin("varfixed",varfixed)
718	write(,)" varfixed = ",varfixed
719
720	write(,)"What is the saturation % of the variable species?"
721	satval=0.8
722	call getin("satval",satval)
723	write(,)" satval = ",satval
724
725
726	! Test of incompatibility:
727	! if varactive, then varfixed should be false
728	if (varactive.and.varfixed) then
729	print*,'if varactive, varfixed must be OFF!'
730	stop
731	endif
732
733	write(,) "Gravitationnal sedimentation ?"
734	sedimentation=.false. ! default value
735	call getin("sedimentation",sedimentation)
736	write(,) " sedimentation = ",sedimentation
737
738	write(,) "Compute water cycle ?"
739	water=.false. ! default value
740	call getin("water",water)
741	write(,) " water = ",water
742
743	! Test of incompatibility:
744	! if water is true, there should be at least a tracer
745	if (water.and.(.not.tracer)) then
746	print*,'if water is ON, tracer must be ON too!'
747	stop
748	endif
749
750	write(,) "Include water condensation ?"
751	watercond=.false. ! default value
752	call getin("watercond",watercond)
753	write(,) " watercond = ",watercond
754
755	! Test of incompatibility:
756	! if watercond is used, then water should be used too
757	if (watercond.and.(.not.water)) then
758	print*,'if watercond is used, water should be used too'
759	stop
760	endif
761
762	write(,) "Include water precipitation ?"
763	waterrain=.false. ! default value
764	call getin("waterrain",waterrain)
765	write(,) " waterrain = ",waterrain
766
767	write(,) "Include surface hydrology ?"
768	hydrology=.false. ! default value
769	call getin("hydrology",hydrology)
770	write(,) " hydrology = ",hydrology
771
772	write(,) "Evolve surface water sources ?"
773	sourceevol=.false. ! default value
774	call getin("sourceevol",sourceevol)
775	write(,) " sourceevol = ",sourceevol
776
777	write(,) "Ice evolution timestep ?"
778	icetstep=100.0 ! default value
779	call getin("icetstep",icetstep)
780	write(,) " icetstep = ",icetstep
781
782	write(,) "Snow albedo ?"
783	albedosnow=0.5 ! default value
784	call getin("albedosnow",albedosnow)
785	write(,) " albedosnow = ",albedosnow
786
787	write(,) "Maximum ice thickness ?"
788	maxicethick=2.0 ! default value
789	call getin("maxicethick",maxicethick)
790	write(,) " maxicethick = ",maxicethick
791
792	write(,) "Freezing point of seawater ?"
793	Tsaldiff=-1.8 ! default value
794	call getin("Tsaldiff",Tsaldiff)
795	write(,) " Tsaldiff = ",Tsaldiff
796
797	write(,) "Does user want to force cpp and mugaz?"
798	force_cpp=.false. ! default value
799	call getin("force_cpp",force_cpp)
800	write(,) " force_cpp = ",force_cpp
801
802	if (force_cpp) then
803	mugaz = -99999.
804	PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
805	call getin("mugaz",mugaz)
806	IF (mugaz.eq.-99999.) THEN
807	PRINT *, "mugaz must be set if force_cpp = T"
808	STOP
809	ELSE
810	write(,) "mugaz=",mugaz
811	ENDIF
812	cpp = -99999.
813	PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
814	call getin("cpp",cpp)
815	IF (cpp.eq.-99999.) THEN
816	PRINT *, "cpp must be set if force_cpp = T"
817	STOP
818	ELSE
819	write(,) "cpp=",cpp
820	ENDIF
821	! else
822	! mugaz=8.314*1000./pr
823	endif
824	call su_gases
825	call calc_cpp_mugaz
826
827	PRINT*,'--------------------------------------------'
828	PRINT*
829	PRINT*
830	ELSE
831	write(,)
832	write(,) 'Cannot read file callphys.def. Is it here ?'
833	stop
834	ENDIF
835
836	8000 FORMAT(t5,a12,l8)
837	8001 FORMAT(t5,a12,i8)
838
839	PRINT*
840	PRINT*,'inifis: daysec',daysec
841	PRINT*
842	PRINT*,'inifis: The radiative transfer is computed:'
843	PRINT*,' each ',iradia,' physical time-step'
844	PRINT,' or each ',iradiadtphys,' seconds'
845	PRINT*
846
847
848	!-----------------------------------------------------------------------
849	! Some more initialization:
850	! ------------------------
851
852	! ALLOCATE ARRAYS IN comgeomfi_h
853	IF (.not. ALLOCATED(lati)) ALLOCATE(lati(ngrid))
854	IF (.not. ALLOCATED(long)) ALLOCATE(long(ngrid))
855	IF (.not. ALLOCATED(area)) ALLOCATE(area(ngrid))
856
857	CALL SCOPY(ngrid,plon,1,long,1)
858	CALL SCOPY(ngrid,plat,1,lati,1)
859	CALL SCOPY(ngrid,parea,1,area,1)
860	totarea=SSUM(ngrid,area,1)
861	call planetwide_sumval(area,totarea_planet)
862
863	!! those are defined in comdiurn_h.F90
864	IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(ngrid))
865	IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(ngrid))
866	IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(ngrid))
867	IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(ngrid))
868
869	DO ig=1,ngrid
870	sinlat(ig)=sin(plat(ig))
871	coslat(ig)=cos(plat(ig))
872	sinlon(ig)=sin(plon(ig))
873	coslon(ig)=cos(plon(ig))
874	ENDDO
875
876	pi=2.*asin(1.) ! NB: pi is a common in comcstfi.h
877
878	! allocate "comsoil_h" arrays
879	call ini_comsoil_h(ngrid)
880
881	END

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

Download in other formats: