Context Navigation

inifis_mod.F90 @ 2219

Last change on this file since 2219 was 2138, checked in by jvatant, 6 years ago
Correct outptuts of dtaui/v with correct ponderations of weights, but now it is attenuation exp(-tau) --JVO
File size: 25.9 KB

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

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

Download in other formats:

Original Format