Context Navigation

inifis.F @ 650

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

Corrected the temperature used to differentiate sublimation and evaporation in watersat_grad
Minor name changes in watercommon
Better physical parametrization of the effective radius of liquid and icy water cloud particles in callcorrk

(for radfixed=true)

Added consistency check in inifis
Moved 1d water initialization from physiqu to rcm1d

File size: 18.4 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
[538]	179	write(,) "Check to see if cpp values used match gases.def ?"
	180	check_cpp_match=.true. ! default value
	181	call getin("check_cpp_match",check_cpp_match)
	182	write(,) " check_cpp_match = ",check_cpp_match
	183
[135]	184	write(,) "call radiative transfer ?"
	185	callrad=.true. ! default value
	186	call getin("callrad",callrad)
	187	write(,) " callrad = ",callrad
	188
	189	write(,) "call correlated-k radiative transfer ?"
	190	corrk=.true. ! default value
	191	call getin("corrk",corrk)
	192	write(,) " corrk = ",corrk
	193
	194	write(,) "call gaseous absorption in the visible bands?",
	195	& "(matters only if callrad=T)"
	196	callgasvis=.false. ! default value
	197	call getin("callgasvis",callgasvis)
	198	write(,) " callgasvis = ",callgasvis
[538]	199
[526]	200	write(,) "call continuum opacities in radiative transfer ?",
	201	& "(matters only if callrad=T)"
[538]	202	Continuum=.true. ! default value
[526]	203	call getin("Continuum",Continuum)
	204	write(,) " Continuum = ",Continuum
[538]	205
	206
[135]	207	write(,) "call turbulent vertical diffusion ?"
	208	calldifv=.true. ! default value
	209	call getin("calldifv",calldifv)
	210	write(,) " calldifv = ",calldifv
	211
[596]	212	write(,) "use turbdiff instead of vdifc ?"
	213	UseTurbDiff=.true. ! default value
	214	call getin("UseTurbDiff",UseTurbDiff)
	215	write(,) " UseTurbDiff = ",UseTurbDiff
	216
[135]	217	write(,) "call convective adjustment ?"
	218	calladj=.true. ! default value
	219	call getin("calladj",calladj)
	220	write(,) " calladj = ",calladj
	221
	222	write(,) "call CO2 condensation ?"
	223	co2cond=.true. ! default value
	224	call getin("co2cond",co2cond)
	225	write(,) " co2cond = ",co2cond
[650]	226	! Test of incompatibility
	227	if (co2cond.and.(.not.tracer)) then
	228	print*,'We need a CO2 ice tracer to condense CO2'
	229	call abort
	230	endif
	231
[253]	232	write(,) "CO2 supersaturation level ?"
	233	co2supsat=1.0 ! default value
	234	call getin("co2supsat",co2supsat)
	235	write(,) " co2supsat = ",co2supsat
	236
	237	write(,) "Radiative timescale for Newtonian cooling ?"
	238	tau_relax=30. ! default value
	239	call getin("tau_relax",tau_relax)
	240	write(,) " tau_relax = ",tau_relax
	241	tau_relax=tau_relax243600 ! convert Earth days --> seconds
	242
[135]	243	write(,)"call thermal conduction in the soil ?"
	244	callsoil=.true. ! default value
	245	call getin("callsoil",callsoil)
	246	write(,) " callsoil = ",callsoil
	247
[253]	248	write(,)"Rad transfer is computed every iradia",
[135]	249	& " physical timestep"
	250	iradia=1 ! default value
	251	call getin("iradia",iradia)
	252	write(,)" iradia = ",iradia
	253
	254	write(,)"Rayleigh scattering ?"
	255	rayleigh=.false.
	256	call getin("rayleigh",rayleigh)
	257	write(,)" rayleigh = ",rayleigh
	258
[253]	259	write(,) "Use blackbody for stellar spectrum ?"
	260	stelbbody=.false. ! default value
	261	call getin("stelbbody",stelbbody)
	262	write(,) " stelbbody = ",stelbbody
[135]	263
[253]	264	write(,) "Stellar blackbody temperature ?"
	265	stelTbb=5800.0 ! default value
	266	call getin("stelTbb",stelTbb)
	267	write(,) " stelTbb = ",stelTbb
	268
[135]	269	write(,)"Output mean OLR in 1D?"
	270	meanOLR=.false.
	271	call getin("meanOLR",meanOLR)
	272	write(,)" meanOLR = ",meanOLR
	273
	274	write(,)"Output spectral OLR in 3D?"
	275	specOLR=.false.
	276	call getin("specOLR",specOLR)
	277	write(,)" specOLR = ",specOLR
	278
[253]	279	write(,)"Operate in kastprof mode?"
	280	kastprof=.false.
	281	call getin("kastprof",kastprof)
	282	write(,)" kastprof = ",kastprof
	283
[596]	284	write(,)"Uniform absorption in radiative transfer?"
[305]	285	graybody=.false.
	286	call getin("graybody",graybody)
	287	write(,)" graybody = ",graybody
[253]	288
	289	! Test of incompatibility:
	290	! if kastprof used, we must be in 1D
	291	if (kastprof.and.(ngridmx.gt.1)) then
	292	print*,'kastprof can only be used in 1D!'
	293	call abort
	294	endif
	295
	296	write(,)"Stratospheric temperature for kastprof mode?"
	297	Tstrat=167.0
	298	call getin("Tstrat",Tstrat)
	299	write(,)" Tstrat = ",Tstrat
	300
	301	write(,)"Remove lower boundary?"
	302	noradsurf=.false.
	303	call getin("noradsurf",noradsurf)
	304	write(,)" noradsurf = ",noradsurf
	305
	306	! Tests of incompatibility:
	307	if (noradsurf.and.callsoil) then
	308	print*,'noradsurf not compatible with soil scheme!'
	309	call abort
	310	endif
[589]	311	!if (noradsurf.and.calldifv) then
	312	! print*,'noradsurf not compatible with a boundary layer!'
	313	! call abort
	314	!endif
[253]	315
	316	write(,)"Use Newtonian cooling for radiative transfer?"
	317	newtonian=.false.
	318	call getin("newtonian",newtonian)
	319	write(,)" newtonian = ",newtonian
	320
	321	! Tests of incompatibility:
	322	if (newtonian.and.corrk) then
	323	print*,'newtonian not compatible with correlated-k!'
	324	call abort
	325	endif
	326	if (newtonian.and.calladj) then
	327	print*,'newtonian not compatible with adjustment!'
	328	call abort
	329	endif
	330	if (newtonian.and.calldifv) then
	331	print*,'newtonian not compatible with a boundary layer!'
	332	call abort
	333	endif
	334
	335	write(,)"Test physics timescale in 1D?"
	336	testradtimes=.false.
	337	call getin("testradtimes",testradtimes)
	338	write(,)" testradtimes = ",testradtimes
	339
	340	! Test of incompatibility:
	341	! if testradtimes used, we must be in 1D
	342	if (testradtimes.and.(ngridmx.gt.1)) then
	343	print*,'testradtimes can only be used in 1D!'
	344	call abort
	345	endif
	346
[135]	347	write(,)"Default planetary temperature?"
	348	tplanet=215.0
	349	call getin("tplanet",tplanet)
	350	write(,)" tplanet = ",tplanet
	351
	352	write(,)"Which star?"
	353	startype=1 ! default value = Sol
	354	call getin("startype",startype)
	355	write(,)" startype = ",startype
	356
	357	write(,)"Value of stellar flux at 1 AU?"
	358	Fat1AU=1356.0 ! default value = Sol today
	359	call getin("Fat1AU",Fat1AU)
	360	write(,)" Fat1AU = ",Fat1AU
	361
	362
	363	! TRACERS:
	364
[253]	365	write(,)"Fixed / inactive aerosol distribution?"
	366	aerofixed=.true. ! default value
[135]	367	call getin("aerofixed",aerofixed)
	368	write(,)" aerofixed = ",aerofixed
	369
[253]	370	write(,)"Varying H2O cloud fraction?"
	371	CLFvarying=.false. ! default value
	372	call getin("CLFvarying",CLFvarying)
	373	write(,)" CLFvarying = ",CLFvarying
	374
	375	write(,)"Value of fixed H2O cloud fraction?"
	376	CLFfixval=1.0 ! default value
	377	call getin("CLFfixval",CLFfixval)
	378	write(,)" CLFfixval = ",CLFfixval
	379
[135]	380	write(,)"Number mixing ratio of CO2 ice particles:"
	381	Nmix_co2=100000. ! default value
	382	call getin("Nmix_co2",Nmix_co2)
	383	write(,)" Nmix_co2 = ",Nmix_co2
	384
	385	write(,)"Number mixing ratio of H2O ice particles:"
	386	Nmix_h2o=10000000. ! default value
	387	call getin("Nmix_h2o",Nmix_h2o)
	388	write(,)" Nmix_h2o = ",Nmix_h2o
	389
[253]	390	write(,)"Opacity of dust (if used):"
	391	dusttau=0. ! default value
	392	call getin("dusttau",dusttau)
	393	write(,)" dusttau = ",dusttau
	394
	395	! Test of incompatibility:
	396	! if less than 3 aerosols, then dusttau should = 0
	397	if ((naerkind.lt.3).and.dusttau.gt.0.) then
	398	print*,'Need naer>2 for radiatively active dust!'
	399	stop
	400	endif
	401
[305]	402	write(,)"Cloud pressure level (with kastprof only):"
	403	cloudlvl=0. ! default value
	404	call getin("cloudlvl",cloudlvl)
	405	write(,)" cloudlvl = ",cloudlvl
	406
[135]	407	write(,)"Is the variable gas species radiatively active?"
[305]	408	Tstrat=167.0
[135]	409	varactive=.false.
	410	call getin("varactive",varactive)
	411	write(,)" varactive = ",varactive
	412
	413	write(,)"Is the variable gas species distribution set?"
	414	varfixed=.false.
	415	call getin("varfixed",varfixed)
	416	write(,)" varfixed = ",varfixed
	417
	418	write(,)"What is the saturation % of the variable species?"
	419	satval=0.8
	420	call getin("satval",satval)
	421	write(,)" satval = ",satval
	422
	423	! Test of incompatibility:
	424	! if no tracers, then aerofixed should be true
	425	if ((.not.tracer).and.(.not.aerofixed)) then
	426	print*,'if tracers are off, aerofixed must be ON!'
	427	stop
	428	endif
	429
	430	! Test of incompatibility:
	431	! if varactive, then varfixed should be false
	432	if (varactive.and.varfixed) then
	433	print*,'if varactive, varfixed must be OFF!'
	434	stop
	435	endif
	436
	437	write(,) "Gravitationnal sedimentation ?"
	438	sedimentation=.true. ! default value
	439	call getin("sedimentation",sedimentation)
	440	write(,) " sedimentation = ",sedimentation
	441
	442	write(,) "Compute water cycle ?"
	443	water=.false. ! default value
	444	call getin("water",water)
	445	write(,) " water = ",water
	446
[622]	447	! Test of incompatibility:
	448	! if water is true, there should be at least a tracer
	449	if (water.and.(.not.tracer)) then
	450	print*,'if water is ON, tracer must be ON too!'
	451	stop
	452	endif
	453
[135]	454	write(,) "Include water condensation ?"
	455	watercond=.false. ! default value
	456	call getin("watercond",watercond)
	457	write(,) " watercond = ",watercond
	458
[622]	459	! Test of incompatibility:
	460	! if watercond is used, then water should be used too
	461	if (watercond.and.(.not.water)) then
	462	print*,'if watercond is used, water should be used too'
	463	stop
	464	endif
	465
[135]	466	write(,) "Include water precipitation ?"
	467	waterrain=.false. ! default value
	468	call getin("waterrain",waterrain)
	469	write(,) " waterrain = ",waterrain
	470
	471	write(,) "Precipitation threshold ?"
	472	rainthreshold=0.011 ! default value (Emmanuel 1997)
	473	call getin("rainthreshold",rainthreshold)
	474	write(,) " rainthreshold = ",rainthreshold
	475
[253]	476	write(,) "Include surface hydrology ?"
	477	hydrology=.false. ! default value
	478	call getin("hydrology",hydrology)
	479	write(,) " hydrology = ",hydrology
[135]	480
[253]	481	write(,) "Evolve surface water sources ?"
	482	sourceevol=.false. ! default value
	483	call getin("sourceevol",sourceevol)
	484	write(,) " sourceevol = ",sourceevol
[135]	485
[486]	486	write(,) "Ice evolution timestep ?"
	487	icetstep=100.0 ! default value
	488	call getin("icetstep",icetstep)
	489	write(,) " icetstep = ",icetstep
	490
[253]	491	write(,) "Snow albedo ?"
	492	albedosnow=0.5 ! default value
	493	call getin("albedosnow",albedosnow)
	494	write(,) " albedosnow = ",albedosnow
[135]	495
[253]	496	write(,) "Maximum ice thickness ?"
	497	maxicethick=2.0 ! default value
	498	call getin("maxicethick",maxicethick)
	499	write(,) " maxicethick = ",maxicethick
[135]	500
[253]	501	write(,) "Freezing point of seawater ?"
	502	Tsaldiff=-1.8 ! default value
	503	call getin("Tsaldiff",Tsaldiff)
	504	write(,) " Tsaldiff = ",Tsaldiff
[135]	505
[589]	506	write(,) "Does user want to force cpp and mugaz?"
	507	force_cpp=.false. ! default value
	508	call getin("force_cpp",force_cpp)
	509	write(,) " force_cpp = ",force_cpp
	510
	511	if (force_cpp) then
	512	mugaz = -99999.
	513	PRINT *,'MEAN MOLECULAR MASS in g mol-1 ?'
	514	call getin("mugaz",mugaz)
	515	IF (mugaz.eq.-99999.) THEN
	516	PRINT *, "mugaz must be set if force_cpp = T"
	517	STOP
	518	ELSE
	519	write(,) "mugaz=",mugaz
	520	ENDIF
	521	cpp = -99999.
	522	PRINT *,'SPECIFIC HEAT CAPACITY in J K-1 kg-1 ?'
	523	call getin("cpp",cpp)
	524	IF (cpp.eq.-99999.) THEN
	525	PRINT *, "cpp must be set if force_cpp = T"
	526	STOP
	527	ELSE
	528	write(,) "cpp=",cpp
	529	ENDIF
	530	else
	531	mugaz=8.314*1000./pr
	532	endif
[253]	533	call su_gases
	534	call calc_cpp_mugaz
[135]	535
	536	PRINT*,'--------------------------------------------'
	537	PRINT*
	538	PRINT*
	539	ELSE
	540	write(,)
	541	write(,) 'Cannot read file callphys.def. Is it here ?'
	542	stop
	543	ENDIF
	544
	545	8000 FORMAT(t5,a12,l8)
	546	8001 FORMAT(t5,a12,i8)
	547
	548	PRINT*
	549	PRINT*,'inifis: daysec',daysec
	550	PRINT*
	551	PRINT*,'inifis: The radiative transfer is computed:'
	552	PRINT*,' each ',iradia,' physical time-step'
	553	PRINT,' or each ',iradiadtphys,' seconds'
	554	PRINT*
	555
	556
	557	!-----------------------------------------------------------------------
	558	! Some more initialization:
	559	! ------------------------
	560
	561	CALL SCOPY(ngrid,plon,1,long,1)
	562	CALL SCOPY(ngrid,plat,1,lati,1)
	563	CALL SCOPY(ngrid,parea,1,area,1)
	564	totarea=SSUM(ngridmx,area,1)
	565
	566	DO ig=1,ngrid
	567	sinlat(ig)=sin(plat(ig))
	568	coslat(ig)=cos(plat(ig))
	569	sinlon(ig)=sin(plon(ig))
	570	coslon(ig)=cos(plon(ig))
	571	ENDDO
	572
	573	pi=2.*asin(1.) ! NB: pi is a common in comcstfi.h
	574
	575	RETURN
	576	END

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

Download in other formats:

Original Format