Context Navigation

inifis.F @ 1217

Last change on this file since 1217 was 1216, checked in by emillour, 11 years ago

Generic model:
Major cleanup, in order to ease the use of LMDZ.GENERIC with (parallel) dynamics
in LMDZ.COMMON: (NB: this will break LMDZ.UNIVERSAL, which should be thrashed
in the near future)

Updated makegcm_* scripts (and makdim) and added the "-full" (to enforce full recomputation of the model) option
In dyn3d: converted control.h to module control_mod.F90 and converted iniadvtrac.F to module infotrac.F90
Added module mod_const_mpi.F90 in dyn3d (not used in serial mode)
Rearanged input/outputs routines everywhere to handle serial/MPI cases. physdem.F => phyredem.F90 , phyetat0.F => phyetat0.F90 ; all read/write routines for startfi files are gathered in module iostart.F90
added parallelism related routines init_phys_lmdz.F90, comgeomphy.F90, dimphy.F90, iniphysiq.F90, mod_grid_phy_lmdz.F90, mod_phys_lmdz_mpi_data.F90, mod_phys_lmdz_mpi_transfert.F90, mod_phys_lmdz_omp_data.F90, mod_phys_lmdz_omp_transfert.F90, mod_phys_lmdz_para.F90, mod_phys_lmdz_transfert_para.F90 in phymars and mod_const_mpi.F90 in dyn3d (for compliance with parallelism)
added created generic routines 'planetwide_maxval' and 'planetwide_minval', in module "planetwide_mod", that enable obtaining the max and min of a field over the whole planet. This should be further imroved with computation of means (possibly area weighed), etc.

File size: 23.3 KB

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

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

Download in other formats:

Original Format