Context Navigation

inifis.F @ 537

Last change on this file since 537 was 526, checked in by jleconte, 13 years ago

13/02/2012 == JL + AS

All outputs are now in netCDF format. Even in 1D (No more G1D)
Clean up of the call to callcorrk when CLFvarying=true
Corrects a bug in writediagspecIR/VI. Output are now in W/m2/cm-1 as a function of the wavenumber in cm-1
Enable writediagspecIR/V to work in the CLFvarying=true case (output now done in Physiq after writediagfi)
Add a simple treatment for the supersaturation of CO2 (see forget et al 2012)
corrects a small bug when no clouds are present in aeropacity

File size: 17.3 KB

Rev	Line
[135]	1	SUBROUTINE inifis(ngrid,nlayer,
	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
[253]	8
[135]	9	!=======================================================================
	10	!
	11	! purpose:
	12	! -------
	13	!
	14	! Initialisation for the physical parametrisations of the LMD
[305]	15	! Generic Model.
[135]	16	!
	17	! author: Frederic Hourdin 15 / 10 /93
	18	! -------
	19	! modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
	20	! Ehouarn Millour (oct. 2008) tracers are now identified
	21	! by their names and may not be contiguously
	22	! stored in the q(:,:,:,:) array
	23	! E.M. (june 2009) use getin routine to load parameters
	24	!
	25	!
	26	! arguments:
	27	! ----------
	28	!
	29	! input:
	30	! ------
	31	!
	32	! ngrid Size of the horizontal grid.
	33	! All internal loops are performed on that grid.
	34	! nlayer Number of vertical layers.
	35	! pdayref Day of reference for the simulation
	36	! pday Number of days counted from the North. Spring
	37	! equinoxe.
	38	!
	39	!=======================================================================
	40	!
	41	!-----------------------------------------------------------------------
	42	! declarations:
	43	! -------------
[374]	44	use datafile_mod, only: datadir
[135]	45	! to use 'getin'
	46	USE ioipsl_getincom
	47	IMPLICIT NONE
	48	#include "dimensions.h"
	49	#include "dimphys.h"
	50	#include "planete.h"
	51	#include "comcstfi.h"
	52	#include "comsaison.h"
	53	#include "comdiurn.h"
	54	#include "comgeomfi.h"
	55	#include "callkeys.h"
	56	#include "surfdat.h"
	57
	58
[305]	59
[135]	60	REAL prad,pg,pr,pcpp,pdaysec,ptimestep
	61
	62	INTEGER ngrid,nlayer
	63	REAL plat(ngrid),plon(ngrid),parea(ngridmx)
	64	integer day_ini
	65	INTEGER ig,ierr
	66
	67	EXTERNAL iniorbit,orbite
	68	EXTERNAL SSUM
	69	REAL SSUM
	70
	71	CHARACTER ch1*12
	72	CHARACTER ch80*80
	73
	74	logical chem, h2o
	75	logical :: parameter, doubleq=.false.
	76
[253]	77	real psurf,pN2 ! added by RW for Gliese 581d N2+CO2
	78
[135]	79	rad=prad
	80	daysec=pdaysec
	81	dtphys=ptimestep
	82	cpp=pcpp
	83	g=pg
	84	r=pr
	85	rcp=r/cpp
	86
	87	avocado = 6.02214179e23 ! added by RW
	88
	89	! --------------------------------------------------------
	90	! The usual Tests
	91	! --------------
	92	IF (nlayer.NE.nlayermx) THEN
	93	PRINT*,'STOP in inifis'
	94	PRINT*,'Probleme de dimensions :'
	95	PRINT*,'nlayer = ',nlayer
	96	PRINT*,'nlayermx = ',nlayermx
	97	STOP
	98	ENDIF
	99
	100	IF (ngrid.NE.ngridmx) THEN
	101	PRINT*,'STOP in inifis'
	102	PRINT*,'Probleme de dimensions :'
	103	PRINT*,'ngrid = ',ngrid
	104	PRINT*,'ngridmx = ',ngridmx
	105	STOP
	106	ENDIF
	107
	108	! --------------------------------------------------------------
	109	! Reading the "callphys.def" file controlling some key options
	110	! --------------------------------------------------------------
	111
	112	! check that 'callphys.def' file is around
	113	OPEN(99,file='callphys.def',status='old',form='formatted'
	114	& ,iostat=ierr)
	115	CLOSE(99)
	116
	117	IF(ierr.EQ.0) THEN
	118	PRINT*
	119	PRINT*
	120	PRINT*,'--------------------------------------------'
	121	PRINT*,' inifis: Parametres pour la physique (callphys.def)'
	122	PRINT*,'--------------------------------------------'
	123
[374]	124	write(,) "Directory where external input files are:"
	125	! default 'datadir' is set in "datadir_mod"
	126	call getin("datadir",datadir) ! default path
	127	write(,) " datadir = ",trim(datadir)
	128
[135]	129	write(,) "Run with or without tracer transport ?"
	130	tracer=.false. ! default value
	131	call getin("tracer",tracer)
	132	write(,) " tracer = ",tracer
	133
	134	write(,) "Diurnal cycle ?"
	135	write(,) "(if diurnal=false, diurnal averaged solar heating)"
	136	diurnal=.true. ! default value
	137	call getin("diurnal",diurnal)
	138	write(,) " diurnal = ",diurnal
	139
	140	write(,) "Seasonal cycle ?"
	141	write(,) "(if season=false, Ls stays constant, to value ",
	142	& "set in 'start'"
	143	season=.true. ! default value
	144	call getin("season",season)
	145	write(,) " season = ",season
	146
	147	write(,) "Tidally resonant rotation ?"
	148	tlocked=.false. ! default value
	149	call getin("tlocked",tlocked)
	150	write(,) "tlocked = ",tlocked
	151
	152	! Test of incompatibility:
	153	! if tlocked, then diurnal should be false
	154	if (tlocked.and.diurnal) then
	155	print*,'If diurnal=true, we should turn off tlocked.'
	156	stop
	157	endif
	158
	159	write(,) "Tidal resonance ratio ?"
	160	nres=0 ! default value
	161	call getin("nres",nres)
	162	write(,) "nres = ",nres
	163
	164	write(,) "Write some extra output to the screen ?"
	165	lwrite=.false. ! default value
	166	call getin("lwrite",lwrite)
	167	write(,) " lwrite = ",lwrite
	168
	169	write(,) "Save statistics in file stats.nc ?"
	170	callstats=.true. ! default value
	171	call getin("callstats",callstats)
	172	write(,) " callstats = ",callstats
	173
	174	write(,) "Test energy conservation of model physics ?"
	175	enertest=.false. ! default value
	176	call getin("enertest",enertest)
	177	write(,) " enertest = ",enertest
	178
	179	write(,) "call radiative transfer ?"
	180	callrad=.true. ! default value
	181	call getin("callrad",callrad)
	182	write(,) " callrad = ",callrad
	183
	184	write(,) "call correlated-k radiative transfer ?"
	185	corrk=.true. ! default value
	186	call getin("corrk",corrk)
	187	write(,) " corrk = ",corrk
	188
	189	write(,) "call gaseous absorption in the visible bands?",
	190	& "(matters only if callrad=T)"
	191	callgasvis=.false. ! default value
	192	call getin("callgasvis",callgasvis)
	193	write(,) " callgasvis = ",callgasvis
[526]	194
	195	write(,) "call continuum opacities in radiative transfer ?",
	196	& "(matters only if callrad=T)"
	197	Continuum=.false. ! default value
	198	call getin("Continuum",Continuum)
	199	write(,) " Continuum = ",Continuum
[135]	200
	201	write(,) "call turbulent vertical diffusion ?"
	202	calldifv=.true. ! default value
	203	call getin("calldifv",calldifv)
	204	write(,) " calldifv = ",calldifv
	205
	206	write(,) "call convective adjustment ?"
	207	calladj=.true. ! default value
	208	call getin("calladj",calladj)
	209	write(,) " calladj = ",calladj
	210
	211	write(,)"Gas is nonideal CO2 ?"
	212	nonideal=.false.
	213	call getin("nonideal",nonideal)
	214	write(,)" nonideal = ",nonideal
	215
	216	! Test of incompatibility
	217	if (enertest.and.nonideal) then
	218	print*,'Energy conservation calculations currently
	219	& assume ideal gas!'
	220	call abort
	221	endif
	222
	223	write(,) "call CO2 condensation ?"
	224	co2cond=.true. ! default value
	225	call getin("co2cond",co2cond)
	226	write(,) " co2cond = ",co2cond
	227
[253]	228	write(,) "CO2 supersaturation level ?"
	229	co2supsat=1.0 ! default value
	230	call getin("co2supsat",co2supsat)
	231	write(,) " co2supsat = ",co2supsat
	232
	233	write(,) "Radiative timescale for Newtonian cooling ?"
	234	tau_relax=30. ! default value
	235	call getin("tau_relax",tau_relax)
	236	write(,) " tau_relax = ",tau_relax
	237	tau_relax=tau_relax243600 ! convert Earth days --> seconds
	238
[135]	239	write(,)"call thermal conduction in the soil ?"
	240	callsoil=.true. ! default value
	241	call getin("callsoil",callsoil)
	242	write(,) " callsoil = ",callsoil
	243
[253]	244	write(,)"Rad transfer is computed every iradia",
[135]	245	& " physical timestep"
	246	iradia=1 ! default value
	247	call getin("iradia",iradia)
	248	write(,)" iradia = ",iradia
	249
	250	write(,)"Rayleigh scattering ?"
	251	rayleigh=.false.
	252	call getin("rayleigh",rayleigh)
	253	write(,)" rayleigh = ",rayleigh
	254
[253]	255	write(,) "Use blackbody for stellar spectrum ?"
	256	stelbbody=.false. ! default value
	257	call getin("stelbbody",stelbbody)
	258	write(,) " stelbbody = ",stelbbody
[135]	259
[253]	260	write(,) "Stellar blackbody temperature ?"
	261	stelTbb=5800.0 ! default value
	262	call getin("stelTbb",stelTbb)
	263	write(,) " stelTbb = ",stelTbb
	264
[135]	265	write(,)"Output mean OLR in 1D?"
	266	meanOLR=.false.
	267	call getin("meanOLR",meanOLR)
	268	write(,)" meanOLR = ",meanOLR
	269
	270	write(,)"Output spectral OLR in 3D?"
	271	specOLR=.false.
	272	call getin("specOLR",specOLR)
	273	write(,)" specOLR = ",specOLR
	274
[253]	275	write(,)"Operate in kastprof mode?"
	276	kastprof=.false.
	277	call getin("kastprof",kastprof)
	278	write(,)" kastprof = ",kastprof
	279
[305]	280	write(,)"Uniform absorption coefficient in IR?"
	281	graybody=.false.
	282	call getin("graybody",graybody)
	283	write(,)" graybody = ",graybody
[253]	284
	285	! Test of incompatibility:
	286	! if kastprof used, we must be in 1D
	287	if (kastprof.and.(ngridmx.gt.1)) then
	288	print*,'kastprof can only be used in 1D!'
	289	call abort
	290	endif
	291
	292	write(,)"Stratospheric temperature for kastprof mode?"
	293	Tstrat=167.0
	294	call getin("Tstrat",Tstrat)
	295	write(,)" Tstrat = ",Tstrat
	296
	297	write(,)"Remove lower boundary?"
	298	noradsurf=.false.
	299	call getin("noradsurf",noradsurf)
	300	write(,)" noradsurf = ",noradsurf
	301
	302	! Tests of incompatibility:
	303	if (noradsurf.and.callsoil) then
	304	print*,'noradsurf not compatible with soil scheme!'
	305	call abort
	306	endif
	307	if (noradsurf.and.calldifv) then
	308	print*,'noradsurf not compatible with a boundary layer!'
	309	call abort
	310	endif
	311
	312
	313	write(,)"Use Newtonian cooling for radiative transfer?"
	314	newtonian=.false.
	315	call getin("newtonian",newtonian)
	316	write(,)" newtonian = ",newtonian
	317
	318	! Tests of incompatibility:
	319	if (newtonian.and.corrk) then
	320	print*,'newtonian not compatible with correlated-k!'
	321	call abort
	322	endif
	323	if (newtonian.and.calladj) then
	324	print*,'newtonian not compatible with adjustment!'
	325	call abort
	326	endif
	327	if (newtonian.and.calldifv) then
	328	print*,'newtonian not compatible with a boundary layer!'
	329	call abort
	330	endif
	331
	332	write(,)"Test physics timescale in 1D?"
	333	testradtimes=.false.
	334	call getin("testradtimes",testradtimes)
	335	write(,)" testradtimes = ",testradtimes
	336
	337	! Test of incompatibility:
	338	! if testradtimes used, we must be in 1D
	339	if (testradtimes.and.(ngridmx.gt.1)) then
	340	print*,'testradtimes can only be used in 1D!'
	341	call abort
	342	endif
	343
[135]	344	write(,)"Default planetary temperature?"
	345	tplanet=215.0
	346	call getin("tplanet",tplanet)
	347	write(,)" tplanet = ",tplanet
	348
	349	write(,)"Which star?"
	350	startype=1 ! default value = Sol
	351	call getin("startype",startype)
	352	write(,)" startype = ",startype
	353
	354	write(,)"Value of stellar flux at 1 AU?"
	355	Fat1AU=1356.0 ! default value = Sol today
	356	call getin("Fat1AU",Fat1AU)
	357	write(,)" Fat1AU = ",Fat1AU
	358
	359
[253]	360	! 1D solar zenith angle
	361	write(,)"Solar zenith angle in 1D?"
	362	szangle=60.0
	363	call getin("szangle",szangle)
	364	write(,)" szangle = ",szangle
	365
[135]	366	! TRACERS:
	367
[253]	368	write(,)"Fixed / inactive aerosol distribution?"
	369	aerofixed=.true. ! default value
[135]	370	call getin("aerofixed",aerofixed)
	371	write(,)" aerofixed = ",aerofixed
	372
[253]	373	write(,)"Varying H2O cloud fraction?"
	374	CLFvarying=.false. ! default value
	375	call getin("CLFvarying",CLFvarying)
	376	write(,)" CLFvarying = ",CLFvarying
	377
	378	write(,)"Value of fixed H2O cloud fraction?"
	379	CLFfixval=1.0 ! default value
	380	call getin("CLFfixval",CLFfixval)
	381	write(,)" CLFfixval = ",CLFfixval
	382
[135]	383	write(,)"Number mixing ratio of CO2 ice particles:"
	384	Nmix_co2=100000. ! default value
	385	call getin("Nmix_co2",Nmix_co2)
	386	write(,)" Nmix_co2 = ",Nmix_co2
	387
	388	write(,)"Number mixing ratio of H2O ice particles:"
	389	Nmix_h2o=10000000. ! default value
	390	call getin("Nmix_h2o",Nmix_h2o)
	391	write(,)" Nmix_h2o = ",Nmix_h2o
	392
[253]	393	write(,)"Opacity of dust (if used):"
	394	dusttau=0. ! default value
	395	call getin("dusttau",dusttau)
	396	write(,)" dusttau = ",dusttau
	397
	398	! Test of incompatibility:
	399	! if less than 3 aerosols, then dusttau should = 0
	400	if ((naerkind.lt.3).and.dusttau.gt.0.) then
	401	print*,'Need naer>2 for radiatively active dust!'
	402	stop
	403	endif
	404
[305]	405	write(,)"Cloud pressure level (with kastprof only):"
	406	cloudlvl=0. ! default value
	407	call getin("cloudlvl",cloudlvl)
	408	write(,)" cloudlvl = ",cloudlvl
	409
[135]	410	write(,)"Is the variable gas species radiatively active?"
[305]	411	Tstrat=167.0
[135]	412	varactive=.false.
	413	call getin("varactive",varactive)
	414	write(,)" varactive = ",varactive
	415
	416	write(,)"Is the variable gas species distribution set?"
	417	varfixed=.false.
	418	call getin("varfixed",varfixed)
	419	write(,)" varfixed = ",varfixed
	420
	421	write(,)"What is the saturation % of the variable species?"
	422	satval=0.8
	423	call getin("satval",satval)
	424	write(,)" satval = ",satval
	425
	426	! Test of incompatibility:
	427	! if no tracers, then aerofixed should be true
	428	if ((.not.tracer).and.(.not.aerofixed)) then
	429	print*,'if tracers are off, aerofixed must be ON!'
	430	stop
	431	endif
	432
	433	! Test of incompatibility:
	434	! if varactive, then varfixed should be false
	435	if (varactive.and.varfixed) then
	436	print*,'if varactive, varfixed must be OFF!'
	437	stop
	438	endif
	439
	440	! Test of incompatibility:
	441
	442	write(,) "Gravitationnal sedimentation ?"
	443	sedimentation=.true. ! default value
	444	call getin("sedimentation",sedimentation)
	445	write(,) " sedimentation = ",sedimentation
	446
	447
	448	! Test of incompatibility:
	449
	450	write(,) "Compute water cycle ?"
	451	water=.false. ! default value
	452	call getin("water",water)
	453	write(,) " water = ",water
	454
	455	write(,) "Include water condensation ?"
	456	watercond=.false. ! default value
	457	call getin("watercond",watercond)
	458	write(,) " watercond = ",watercond
	459
	460	write(,) "Include water precipitation ?"
	461	waterrain=.false. ! default value
	462	call getin("waterrain",waterrain)
	463	write(,) " waterrain = ",waterrain
	464
	465	write(,) "Precipitation threshold ?"
	466	rainthreshold=0.011 ! default value (Emmanuel 1997)
	467	call getin("rainthreshold",rainthreshold)
	468	write(,) " rainthreshold = ",rainthreshold
	469
[253]	470	write(,) "Include surface hydrology ?"
	471	hydrology=.false. ! default value
	472	call getin("hydrology",hydrology)
	473	write(,) " hydrology = ",hydrology
[135]	474
[253]	475	write(,) "Evolve surface water sources ?"
	476	sourceevol=.false. ! default value
	477	call getin("sourceevol",sourceevol)
	478	write(,) " sourceevol = ",sourceevol
[135]	479
[486]	480	write(,) "Ice evolution timestep ?"
	481	icetstep=100.0 ! default value
	482	call getin("icetstep",icetstep)
	483	write(,) " icetstep = ",icetstep
	484
[253]	485	write(,) "Snow albedo ?"
	486	albedosnow=0.5 ! default value
	487	call getin("albedosnow",albedosnow)
	488	write(,) " albedosnow = ",albedosnow
[135]	489
[253]	490	write(,) "Maximum ice thickness ?"
	491	maxicethick=2.0 ! default value
	492	call getin("maxicethick",maxicethick)
	493	write(,) " maxicethick = ",maxicethick
[135]	494
[253]	495	write(,) "Freezing point of seawater ?"
	496	Tsaldiff=-1.8 ! default value
	497	call getin("Tsaldiff",Tsaldiff)
	498	write(,) " Tsaldiff = ",Tsaldiff
[135]	499
	500	! Test of incompatibility:
[253]	501	! if watercond is used, then water should be used too
[135]	502
[253]	503	if (watercond.and.(.not.watercond)) then
	504	print*,'if watercond is used, water should be used too'
[135]	505	stop
	506	endif
	507
[253]	508	mugaz=0.
	509	cpp=0.
	510	call su_gases
	511	call calc_cpp_mugaz
[135]	512
	513	PRINT*,'--------------------------------------------'
	514	PRINT*
	515	PRINT*
	516	ELSE
	517	write(,)
	518	write(,) 'Cannot read file callphys.def. Is it here ?'
	519	stop
	520	ENDIF
	521
	522	8000 FORMAT(t5,a12,l8)
	523	8001 FORMAT(t5,a12,i8)
	524
	525	PRINT*
	526	PRINT*,'inifis: daysec',daysec
	527	PRINT*
	528	PRINT*,'inifis: The radiative transfer is computed:'
	529	PRINT*,' each ',iradia,' physical time-step'
	530	PRINT,' or each ',iradiadtphys,' seconds'
	531	PRINT*
	532
	533
	534	!-----------------------------------------------------------------------
	535	! Some more initialization:
	536	! ------------------------
	537
	538	CALL SCOPY(ngrid,plon,1,long,1)
	539	CALL SCOPY(ngrid,plat,1,lati,1)
	540	CALL SCOPY(ngrid,parea,1,area,1)
	541	totarea=SSUM(ngridmx,area,1)
	542
	543	DO ig=1,ngrid
	544	sinlat(ig)=sin(plat(ig))
	545	coslat(ig)=cos(plat(ig))
	546	sinlon(ig)=sin(plon(ig))
	547	coslon(ig)=cos(plon(ig))
	548	ENDDO
	549
	550	pi=2.*asin(1.) ! NB: pi is a common in comcstfi.h
	551
	552	RETURN
	553	END

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

Download in other formats:

Original Format