Context Navigation

inifis.F @ 1421

Last change on this file since 1421 was 1419, checked in by jleconte, 10 years ago
Correct radiative transfer bug for the Top layer in callcorrk following commit by SG.
File size: 24.4 KB

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

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

Download in other formats:

Original Format