Context Navigation

source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/inifis.F @ 298

Last change on this file since 298 was 298, checked in by milmd, 10 years ago
Less output messages are written. On 20000 cores it is better. In LMDZ, only master of MPI and OpenMP can write.
File size: 28.0 KB

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

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

Download in other formats: