Context Navigation

inifis.F @ 937

Last change on this file since 937 was 873, checked in by aslmd, 12 years ago

LMDZ.GENERIC (and LMDZ.UNIVERSAL)

Optimized calculations for continuum (done for H2 and He, to be done for others)
- new common bilinear interpolation routine (bilinearbig)
- optimization: only one calculation is actually needed

to find indexes of wavelength for bilinear interpolation
... because this will not change with level and integration step!

optimization: use while loop in bilinearbig
completely similar results obtained (test case for a gas giant, many simulated days)

NB: those changes really improve gcm speed (factor 2.2 for whole model!)

continuum was very expensive, now very cheap
--> e.g. 1 day, 25 dyn ts, 5 phys ts
--> before: 243 seconds (including 120 seconds for continuum bilinear interpolation)
--> after: 108 seconds

Corrected a bug: Continuum in inifis instead of continuum

... until now, most users (unbeknownst to them) were running with the continuum by default!

Cosmetic changes in optcv (mostly spaces and line breaks)

... so that comparisons with optci are easy e.g. through vimdiff

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

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

Download in other formats:

Original Format