Context Navigation

inifis.F @ 545

Last change on this file since 545 was 544, checked in by acolaitis, 13 years ago

27/02/12 == AC

Continuation of thermals setting, comparisons with mesoscale results on Case C
Added possibility to call gcm (or 1d) with constant prescribed sensible heat flux, in the spirit of direct comparisons with LES

... This is directly comparable to the variable tke_heat_flux in namelist.input
... Requires the use of yamada4.F from terrestrial GCM (mixes more, seem more numerically stable)
... Usually requires high timesteps (>1000) to avoids crashes. Best approach is to compare
height of first model level z1 and teta_1 between LES and 1D, and increase the timesteps until results
between the two models are comparable (might require a slitghly different tke_heat_flux between the two models
due to difference in vertical diffusion schemes and subgrid effects)
... Syntax for use is to add "tke_heat_flux = ..." in callphys.def

Corrected some stuff with tke transport in thermals (which is anyway deactivated but one day maybe...)

File size: 26.1 KB

Rev	Line
[234]	1	SUBROUTINE inifis(
[226]	2	$ ngrid,nlayer
	3	$ ,day_ini,pdaysec,ptimestep
	4	$ ,plat,plon,parea
	5	$ ,prad,pg,pr,pcpp
	6	#ifdef MESOSCALE
	7	#include "meso_inc/meso_inc_inifisinvar.F"
	8	#endif
	9	$ )
[42]	10	!
	11	!=======================================================================
	12	!
	13	! purpose:
	14	! -------
	15	!
	16	! Initialisation for the physical parametrisations of the LMD
	17	! martian atmospheric general circulation modele.
	18	!
	19	! author: Frederic Hourdin 15 / 10 /93
	20	! -------
	21	! modified: Sebastien Lebonnois 11/06/2003 (new callphys.def)
	22	! Ehouarn Millour (oct. 2008) tracers are now identified
	23	! by their names and may not be contiguously
	24	! stored in the q(:,:,:,:) array
	25	! E.M. (june 2009) use getin routine to load parameters
[234]	26	! adapted to the mesoscale use - Aymeric Spiga - 01/2007-07/2011
[42]	27	!
	28	!
	29	! arguments:
	30	! ----------
	31	!
	32	! input:
	33	! ------
	34	!
	35	! ngrid Size of the horizontal grid.
	36	! All internal loops are performed on that grid.
	37	! nlayer Number of vertical layers.
	38	! pdayref Day of reference for the simulation
	39	! pday Number of days counted from the North. Spring
	40	! equinoxe.
	41	!
	42	!=======================================================================
	43	!
	44	!-----------------------------------------------------------------------
	45	! declarations:
	46	! -------------
	47	! to use 'getin'
	48	USE ioipsl_getincom
	49	IMPLICIT NONE
	50	#include "dimensions.h"
	51	#include "dimphys.h"
	52	#include "planete.h"
	53	#include "comcstfi.h"
	54	#include "comsaison.h"
	55	#include "comdiurn.h"
	56	#include "comgeomfi.h"
	57	#include "callkeys.h"
	58	#include "surfdat.h"
[226]	59	#include "dimradmars.h"
	60	#include "yomaer.h"
[42]	61	#include "datafile.h"
[234]	62	#include "slope.h"
[358]	63	#include "microphys.h"
[420]	64	#include "tracer.h"
[226]	65	#ifdef MESOSCALE
	66	#include "comsoil.h" !!MESOSCALE -- needed to fill volcapa
	67	#include "meso_inc/meso_inc_inifisvar.F"
	68	#endif
[42]	69	REAL prad,pg,pr,pcpp,pdaysec
[234]	70
[226]	71	REAL ptimestep
	72	INTEGER day_ini
[234]	73
[226]	74	INTEGER ngrid,nlayer
[42]	75	REAL plat(ngrid),plon(ngrid),parea(ngridmx)
	76	INTEGER ig,ierr
	77
	78	! EXTERNAL iniorbit,orbite
	79	EXTERNAL SSUM
	80	REAL SSUM
	81
	82	CHARACTER ch1*12
	83	CHARACTER ch80*80
	84
	85	! logical chem, h2o
	86
	87	! chem = .false.
	88	! h2o = .false.
	89
[226]	90	rad=prad
[42]	91	cpp=pcpp
	92	g=pg
	93	r=pr
	94	rcp=r/cpp
[226]	95	daysec=pdaysec
	96	dtphys=ptimestep
[234]	97	#ifdef MESOSCALE
[226]	98	#include "meso_inc/meso_inc_inifisini.F"
	99	#endif
[42]	100
	101	! --------------------------------------------------------
	102	! The usual Tests
	103	! --------------
	104	IF (nlayer.NE.nlayermx) THEN
	105	PRINT*,'STOP in inifis'
	106	PRINT*,'Probleme de dimensions :'
	107	PRINT*,'nlayer = ',nlayer
	108	PRINT*,'nlayermx = ',nlayermx
	109	STOP
	110	ENDIF
	111
	112	IF (ngrid.NE.ngridmx) THEN
	113	PRINT*,'STOP in inifis'
	114	PRINT*,'Probleme de dimensions :'
	115	PRINT*,'ngrid = ',ngrid
	116	PRINT*,'ngridmx = ',ngridmx
	117	STOP
	118	ENDIF
	119
	120	! --------------------------------------------------------------
	121	! Reading the "callphys.def" file controlling some key options
	122	! --------------------------------------------------------------
	123
	124	! check that 'callphys.def' file is around
	125	OPEN(99,file='callphys.def',status='old',form='formatted'
	126	& ,iostat=ierr)
	127	CLOSE(99)
	128
	129	IF(ierr.EQ.0) THEN
	130	PRINT*
	131	PRINT*
	132	PRINT*,'--------------------------------------------'
	133	PRINT*,' inifis: Parameters for the physics (callphys.def)'
	134	PRINT*,'--------------------------------------------'
	135
	136	write(,) "Directory where external input files are:"
[426]	137	datafile="/u/forget/WWW/datagcm/datafile"
[226]	138	call getin("datadir",datafile) ! default path
[42]	139	write(,) " datafile = ",trim(datafile)
	140
	141	write(,) "Run with or without tracer transport ?"
	142	tracer=.false. ! default value
	143	call getin("tracer",tracer)
	144	write(,) " tracer = ",tracer
	145
	146	write(,) "Diurnal cycle ?"
	147	write(,) "(if diurnal=False, diurnal averaged solar heating)"
	148	diurnal=.true. ! default value
	149	call getin("diurnal",diurnal)
	150	write(,) " diurnal = ",diurnal
	151
	152	write(,) "Seasonal cycle ?"
	153	write(,) "(if season=False, Ls stays constant, to value ",
	154	& "set in 'start'"
	155	season=.true. ! default value
	156	call getin("season",season)
	157	write(,) " season = ",season
	158
	159	write(,) "Write some extra output to the screen ?"
	160	lwrite=.false. ! default value
	161	call getin("lwrite",lwrite)
	162	write(,) " lwrite = ",lwrite
	163
	164	write(,) "Save statistics in file stats.nc ?"
[226]	165	#ifdef MESOSCALE
[42]	166	callstats=.false. ! default value
[226]	167	#else
	168	callstats=.true. ! default value
	169	#endif
[42]	170	call getin("callstats",callstats)
	171	write(,) " callstats = ",callstats
	172
	173	write(,) "Save EOF profiles in file 'profiles' for ",
	174	& "Climate Database?"
	175	calleofdump=.false. ! default value
	176	call getin("calleofdump",calleofdump)
	177	write(,) " calleofdump = ",calleofdump
	178
	179	write(,) "Dust scenario: 1=constant dust (read from startfi",
	180	& " or set as tauvis); 2=Viking scenario; =3 MGS scenario,",
	181	& "=4 Mars Year 24 from TES assimilation, ",
	182	& "=24,25 or 26 :Mars Year 24,25 or 26 from TES assimilation"
	183	iaervar=3 ! default value
	184	call getin("iaervar",iaervar)
	185	write(,) " iaervar = ",iaervar
	186
	187	write(,) "Reference (visible) dust opacity at 700 Pa ",
	188	& "(matters only if iaervar=1)"
	189	! NB: default value of tauvis is set/read in startfi.nc file
	190	call getin("tauvis",tauvis)
	191	write(,) " tauvis = ",tauvis
	192
	193	write(,) "Dust vertical distribution:"
	194	write(,) "(=1 top set by topdustref parameter;",
	195	& " =2 Viking scenario; =3 MGS scenario)"
	196	iddist=3 ! default value
	197	call getin("iddist",iddist)
	198	write(,) " iddist = ",iddist
	199
	200	write(,) "Dust top altitude (km). (Matters only if iddist=1)"
	201	topdustref= 90.0 ! default value
	202	call getin("topdustref",topdustref)
	203	write(,) " topdustref = ",topdustref
	204
[544]	205	write(,) "Prescribed surface thermal flux (H/(rho*cp),K m/s)"
	206	tke_heat_flux=0. ! default value
	207	call getin("tke_heat_flux",tke_heat_flux)
	208	write(,) " tke_heat_flux = ",tke_heat_flux
	209	write(,) " 0 means the usual schemes are computing"
	210
[42]	211	write(,) "call radiative transfer ?"
	212	callrad=.true. ! default value
	213	call getin("callrad",callrad)
	214	write(,) " callrad = ",callrad
	215
[234]	216	write(,) "call slope insolation scheme ?",
	217	& "(matters only if callrad=T)"
	218	#ifdef MESOSCALE
	219	callslope=.true. ! default value
	220	#else
	221	callslope=.false. ! default value (not supported yet)
	222	#endif
	223	call getin("callslope",callslope)
	224	write(,) " callslope = ",callslope
	225
[42]	226	write(,) "call NLTE radiative schemes ?",
	227	& "(matters only if callrad=T)"
	228	callnlte=.false. ! default value
	229	call getin("callnlte",callnlte)
	230	write(,) " callnlte = ",callnlte
	231
[414]	232	nltemodel=0 !default value
	233	write(,) "NLTE model?"
	234	write(,) "0 -> old model, static O"
	235	write(,) "1 -> old model, dynamic O"
	236	write(,) "2 -> new model"
	237	write(,) "(matters only if callnlte=T)"
	238	call getin("nltemodel",nltemodel)
	239	write(,) " nltemodel = ",nltemodel
	240
[42]	241	write(,) "call CO2 NIR absorption ?",
	242	& "(matters only if callrad=T)"
	243	callnirco2=.false. ! default value
	244	call getin("callnirco2",callnirco2)
	245	write(,) " callnirco2 = ",callnirco2
[414]	246
	247	write(,) "New NIR NLTE correction ?",
	248	$ "0-> old model (no correction)",
	249	$ "1-> new correction",
	250	$ "(matters only if callnirco2=T)"
	251	call getin("nircorr",nircorr)
	252	write(,) " nircorr = ",nircorr
	253
[42]	254	write(,) "call turbulent vertical diffusion ?"
	255	calldifv=.true. ! default value
	256	call getin("calldifv",calldifv)
	257	write(,) " calldifv = ",calldifv
	258
[162]	259	write(,) "call thermals ?"
	260	calltherm=.false. ! default value
	261	call getin("calltherm",calltherm)
	262	write(,) " calltherm = ",calltherm
	263
	264	write(,) "output thermal diagnostics ?"
	265	outptherm=.false. ! default value
	266	call getin("outptherm",outptherm)
	267	write(,) " outptherm = ",outptherm
	268
[42]	269	write(,) "call convective adjustment ?"
	270	calladj=.true. ! default value
	271	call getin("calladj",calladj)
	272	write(,) " calladj = ",calladj
	273
[162]	274	if (calltherm .and. (.not. calladj)) then
	275	print*,'Convadj has to be activated when using thermals'
	276	stop
	277	endif
[42]	278
[284]	279	write(,) "call Richardson-based surface layer ?"
	280	callrichsl=.false. ! default value
	281	call getin("callrichsl",callrichsl)
	282	write(,) " callrichsl = ",callrichsl
	283
[288]	284	if (calltherm .and. .not.callrichsl) then
	285	print*,'WARNING WARNING WARNING'
	286	print*,'if calltherm=T we strongly advise that '
	287	print*,'you use the new surface layer scheme '
	288	print*,'by setting callrichsl=T '
	289	endif
	290
[42]	291	write(,) "call CO2 condensation ?"
	292	callcond=.true. ! default value
	293	call getin("callcond",callcond)
	294	write(,) " callcond = ",callcond
	295
	296	write(,)"call thermal conduction in the soil ?"
	297	callsoil=.true. ! default value
	298	call getin("callsoil",callsoil)
	299	write(,) " callsoil = ",callsoil
	300
	301
	302	write(,)"call Lott's gravity wave/subgrid topography ",
	303	& "scheme ?"
	304	calllott=.true. ! default value
	305	call getin("calllott",calllott)
	306	write(,)" calllott = ",calllott
	307
	308
	309	write(,)"rad.transfer is computed every iradia",
	310	& " physical timestep"
	311	iradia=1 ! default value
	312	call getin("iradia",iradia)
	313	write(,)" iradia = ",iradia
	314
	315
	316	write(,)"Output of the exchange coefficient mattrix ?",
	317	& "(for diagnostics only)"
	318	callg2d=.false. ! default value
	319	call getin("callg2d",callg2d)
	320	write(,)" callg2d = ",callg2d
	321
	322	write(,)"Rayleigh scattering : (should be .false. for now)"
	323	rayleigh=.false.
	324	call getin("rayleigh",rayleigh)
	325	write(,)" rayleigh = ",rayleigh
	326
	327
	328	! TRACERS:
	329
[358]	330	! dustbin
[42]	331	write(,)"Transported dust ? (if >0, use 'dustbin' dust bins)"
	332	dustbin=0 ! default value
	333	call getin("dustbin",dustbin)
	334	write(,)" dustbin = ",dustbin
[358]	335	! active
[42]	336	write(,)"Radiatively active dust ? (matters if dustbin>0)"
	337	active=.false. ! default value
	338	call getin("active",active)
	339	write(,)" active = ",active
	340
	341	! Test of incompatibility:
	342	! if active is used, then dustbin should be > 0
	343
	344	if (active.and.(dustbin.lt.1)) then
	345	print*,'if active is used, then dustbin should > 0'
	346	stop
	347	endif
[358]	348	! doubleq
[42]	349	write(,)"use mass and number mixing ratios to predict",
	350	& " dust size ?"
	351	doubleq=.false. ! default value
	352	call getin("doubleq",doubleq)
	353	write(,)" doubleq = ",doubleq
[358]	354	! submicron
[42]	355	submicron=.false. ! default value
	356	call getin("submicron",submicron)
	357	write(,)" submicron = ",submicron
	358
	359	! Test of incompatibility:
	360	! if doubleq is used, then dustbin should be 2
	361
	362	if (doubleq.and.(dustbin.ne.2)) then
	363	print*,'if doubleq is used, then dustbin should be 2'
	364	stop
	365	endif
	366	if (doubleq.and.submicron.and.(nqmx.LT.3)) then
	367	print*,'If doubleq is used with a submicron tracer,'
	368	print*,' then the number of tracers has to be'
	369	print*,' larger than 3.'
	370	stop
	371	endif
[358]	372	! lifting
[42]	373	write(,)"dust lifted by GCM surface winds ?"
	374	lifting=.false. ! default value
	375	call getin("lifting",lifting)
	376	write(,)" lifting = ",lifting
	377
	378	! Test of incompatibility:
	379	! if lifting is used, then dustbin should be > 0
	380
	381	if (lifting.and.(dustbin.lt.1)) then
	382	print*,'if lifting is used, then dustbin should > 0'
	383	stop
	384	endif
[358]	385	! callddevil
[42]	386	write(,)" dust lifted by dust devils ?"
	387	callddevil=.false. !default value
	388	call getin("callddevil",callddevil)
	389	write(,)" callddevil = ",callddevil
	390
	391	! Test of incompatibility:
	392	! if dustdevil is used, then dustbin should be > 0
	393
	394	if (callddevil.and.(dustbin.lt.1)) then
	395	print*,'if dustdevil is used, then dustbin should > 0'
	396	stop
	397	endif
[358]	398	! sedimentation
[42]	399	write(,) "Gravitationnal sedimentation ?"
	400	sedimentation=.true. ! default value
	401	call getin("sedimentation",sedimentation)
	402	write(,) " sedimentation = ",sedimentation
[358]	403	! activice
[42]	404	write(,) "Radiatively active transported atmospheric ",
	405	& "water ice ?"
	406	activice=.false. ! default value
	407	call getin("activice",activice)
	408	write(,) " activice = ",activice
[358]	409	! water
[42]	410	write(,) "Compute water cycle ?"
	411	water=.false. ! default value
	412	call getin("water",water)
	413	write(,) " water = ",water
	414
	415	! Test of incompatibility:
	416
	417	if (activice.and..not.water) then
	418	print*,'if activice is used, water should be used too'
	419	stop
	420	endif
	421
	422	if (water.and..not.tracer) then
	423	print*,'if water is used, tracer should be used too'
	424	stop
	425	endif
[420]	426
[455]	427	! water ice clouds effective variance distribution for sedimentaion
	428	write(,) "effective variance for water ice clouds ?"
	429	nuice_sed=0.45
	430	call getin("nuice_sed",nuice_sed)
	431	write(,) "water_param nueff Sedimentation:", nuice_sed
	432
[420]	433	! ccn factor if no scavenging
	434	write(,) "water param CCN reduc. factor ?", ccn_factor
	435	ccn_factor = 4.5
	436	call getin("ccn_factor",ccn_factor)
	437	write(,)" ccn_factor = ",ccn_factor
	438	write(,)"Careful: only used when microphys=F, otherwise"
	439	write(,)"the contact parameter is used instead;"
	440
[358]	441	! microphys
	442	write(,)"Microphysical scheme for water-ice clouds?"
	443	microphys=.false. ! default value
	444	call getin("microphys",microphys)
	445	write(,)" microphys = ",microphys
[42]	446
[358]	447	! microphysical parameter contact
	448	write(,) "water contact parameter ?"
	449	mteta = 0.95
	450	call getin("mteta",mteta)
	451	write(,) "mteta = ", mteta
	452
	453	! scavenging
	454	write(,)"Dust scavenging by H2O/CO2 snowfall ?"
	455	scavenging=.false. ! default value
	456	call getin("scavenging",scavenging)
	457	write(,)" scavenging = ",scavenging
	458
	459
[42]	460	! Test of incompatibility:
[358]	461	! if scavenging is used, then dustbin should be > 0
[42]	462
[358]	463	if (scavenging.and.(dustbin.lt.1)) then
	464	print*,'if scavenging is used, then dustbin should > 0'
	465	stop
	466	endif
	467	if ((microphys.and..not.doubleq).or.
	468	& (microphys.and..not.active).or.
	469	& (microphys.and..not.scavenging).or.
	470	& (microphys.and..not.water)) then
	471	print*,'if microphys is used, then doubleq,'
	472	print*,'active, water and scavenging must all be used!'
	473	stop
	474	endif
	475	if ((scavenging.and..not.doubleq).or.
	476	& (scavenging.and..not.active).or.
	477	& (scavenging.and..not.water).or.
	478	& (scavenging.and..not.microphys)) then
	479	print*,'if scavenging is used, then doubleq,'
	480	print*,'active, water and microphys must all be used!'
	481	stop
	482	endif
	483
	484	! Test of incompatibility:
	485
[42]	486	write(,) "Permanent water caps at poles ?",
	487	& " .true. is RECOMMENDED"
	488	write(,) "(with .true., North cap is a source of water ",
	489	& "and South pole is a cold trap)"
	490	caps=.true. ! default value
	491	call getin("caps",caps)
	492	write(,) " caps = ",caps
	493
[358]	494	! albedo_h2o_ice
[283]	495	write(,) "water ice albedo ?"
	496	albedo_h2o_ice=0.45
	497	call getin("albedo_h2o_ice",albedo_h2o_ice)
	498	write(,) " albedo_h2o_ice = ",albedo_h2o_ice
[358]	499	! inert_h2o_ice
[283]	500	write(,) "water ice thermal inertia ?"
	501	inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
	502	call getin("inert_h2o_ice",inert_h2o_ice)
	503	write(,) " inert_h2o_ice = ",inert_h2o_ice
[358]	504	! frost_albedo_threshold
[283]	505	write(,) "frost thickness threshold for albedo ?"
[285]	506	frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2)
[283]	507	call getin("frost_albedo_threshold",
	508	& frost_albedo_threshold)
	509	write(,) " frost_albedo_threshold = ",
	510	& frost_albedo_threshold
	511
[485]	512	! call Titus crocus line -- DEFAULT IS NONE
	513	write(,) "Titus crocus line ?"
	514	tituscap=.false. ! default value
	515	call getin("tituscap",tituscap)
	516	write(,) "tituscap",tituscap
	517
[283]	518	!!!!!!!!!!!!!!!! TEMP !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	519	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	520	write(,) "temp tag for water caps ?"
	521	temptag=.false.
	522	call getin("temptag",temptag)
	523	write(,) " temptag = ",temptag
	524	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	525	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	526
[42]	527	write(,) "photochemistry: include chemical species"
	528	photochem=.false. ! default value
	529	call getin("photochem",photochem)
	530	write(,) " photochem = ",photochem
	531
	532
	533	! THERMOSPHERE
	534
	535	write(,) "call thermosphere ?"
	536	callthermos=.false. ! default value
	537	call getin("callthermos",callthermos)
	538	write(,) " callthermos = ",callthermos
	539
	540
	541	write(,) " water included without cycle ",
	542	& "(only if water=.false.)"
	543	thermoswater=.false. ! default value
	544	call getin("thermoswater",thermoswater)
	545	write(,) " thermoswater = ",thermoswater
	546
	547	write(,) "call thermal conduction ?",
	548	& " (only if callthermos=.true.)"
	549	callconduct=.false. ! default value
	550	call getin("callconduct",callconduct)
	551	write(,) " callconduct = ",callconduct
	552
	553	write(,) "call EUV heating ?",
	554	& " (only if callthermos=.true.)"
	555	calleuv=.false. ! default value
	556	call getin("calleuv",calleuv)
	557	write(,) " calleuv = ",calleuv
	558
	559	write(,) "call molecular viscosity ?",
	560	& " (only if callthermos=.true.)"
	561	callmolvis=.false. ! default value
	562	call getin("callmolvis",callmolvis)
	563	write(,) " callmolvis = ",callmolvis
	564
	565	write(,) "call molecular diffusion ?",
	566	& " (only if callthermos=.true.)"
	567	callmoldiff=.false. ! default value
	568	call getin("callmoldiff",callmoldiff)
	569	write(,) " callmoldiff = ",callmoldiff
	570
	571
	572	write(,) "call thermospheric photochemistry ?",
	573	& " (only if callthermos=.true.)"
	574	thermochem=.false. ! default value
	575	call getin("thermochem",thermochem)
	576	write(,) " thermochem = ",thermochem
	577
	578	write(,) "date for solar flux calculation:",
	579	& " (1985 < date < 2002)"
	580	write(,) "(Solar min=1996.4 ave=1993.4 max=1990.6)"
	581	solarcondate=1993.4 ! default value
	582	call getin("solarcondate",solarcondate)
	583	write(,) " solarcondate = ",solarcondate
	584
	585
	586	if (.not.callthermos) then
	587	if (thermoswater) then
	588	print*,'if thermoswater is set, callthermos must be true'
	589	stop
	590	endif
	591	if (callconduct) then
	592	print*,'if callconduct is set, callthermos must be true'
	593	stop
	594	endif
	595	if (calleuv) then
	596	print*,'if calleuv is set, callthermos must be true'
	597	stop
	598	endif
	599	if (callmolvis) then
	600	print*,'if callmolvis is set, callthermos must be true'
	601	stop
	602	endif
	603	if (callmoldiff) then
	604	print*,'if callmoldiff is set, callthermos must be true'
	605	stop
	606	endif
	607	if (thermochem) then
	608	print*,'if thermochem is set, callthermos must be true'
	609	stop
	610	endif
	611	endif
	612
	613	! Test of incompatibility:
	614	! if photochem is used, then water should be used too
	615
	616	if (photochem.and..not.water) then
	617	print*,'if photochem is used, water should be used too'
	618	stop
	619	endif
	620
	621	! if callthermos is used, then thermoswater should be used too
	622	! (if water not used already)
	623
	624	if (callthermos .and. .not.water) then
	625	if (callthermos .and. .not.thermoswater) then
	626	print*,'if callthermos is used, water or thermoswater
	627	& should be used too'
	628	stop
	629	endif
	630	endif
	631
	632	PRINT*,'--------------------------------------------'
	633	PRINT*
	634	PRINT*
	635	ELSE
	636	write(,)
	637	write(,) 'Cannot read file callphys.def. Is it here ?'
	638	stop
	639	ENDIF
	640
	641	8000 FORMAT(t5,a12,l8)
	642	8001 FORMAT(t5,a12,i8)
	643
	644	PRINT*
	645	PRINT*,'inifis: daysec',daysec
	646	PRINT*
	647	PRINT*,'inifis: The radiative transfer is computed:'
	648	PRINT*,' each ',iradia,' physical time-step'
	649	PRINT,' or each ',iradiadtphys,' seconds'
	650	PRINT*
	651	! --------------------------------------------------------------
	652	! Managing the Longwave radiative transfer
	653	! --------------------------------------------------------------
	654
	655	! In most cases, the run just use the following values :
	656	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	657	callemis=.true.
	658	! ilwd=10*int(daysec/dtphys) ! bug before 22/10/01
	659	ilwd=1
	660	ilwn=1 !2
	661	ilwb=1 !2
	662	linear=.true.
	663	ncouche=3
	664	alphan=0.4
	665	semi=0
	666
	667	! BUT people working hard on the LW may want to read them in 'radia.def'
	668	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	669
	670	OPEN(99,file='radia.def',status='old',form='formatted'
	671	. ,iostat=ierr)
	672	IF(ierr.EQ.0) THEN
	673	write(,) 'inifis: Reading radia.def !!!'
	674	READ(99,fmt='(a)') ch1
	675	READ(99,*) callemis
	676	WRITE(*,8000) ch1,callemis
	677
	678	READ(99,fmt='(a)') ch1
	679	READ(99,*) iradia
	680	WRITE(*,8001) ch1,iradia
	681
	682	READ(99,fmt='(a)') ch1
	683	READ(99,*) ilwd
	684	WRITE(*,8001) ch1,ilwd
	685
	686	READ(99,fmt='(a)') ch1
	687	READ(99,*) ilwn
	688	WRITE(*,8001) ch1,ilwn
	689
	690	READ(99,fmt='(a)') ch1
	691	READ(99,*) linear
	692	WRITE(*,8000) ch1,linear
	693
	694	READ(99,fmt='(a)') ch1
	695	READ(99,*) ncouche
	696	WRITE(*,8001) ch1,ncouche
	697
	698	READ(99,fmt='(a)') ch1
	699	READ(99,*) alphan
	700	WRITE(,) ch1,alphan
	701
	702	READ(99,fmt='(a)') ch1
	703	READ(99,*) ilwb
	704	WRITE(*,8001) ch1,ilwb
	705
	706
	707	READ(99,fmt='(a)') ch1
	708	READ(99,'(l1)') callg2d
	709	WRITE(*,8000) ch1,callg2d
	710
	711	READ(99,fmt='(a)') ch1
	712	READ(99,*) semi
	713	WRITE(,) ch1,semi
	714	end if
	715	CLOSE(99)
	716
	717	!-----------------------------------------------------------------------
	718	! Some more initialization:
	719	! ------------------------
	720
[226]	721	! in 'comgeomfi.h'
[42]	722	CALL SCOPY(ngrid,plon,1,long,1)
	723	CALL SCOPY(ngrid,plat,1,lati,1)
	724	CALL SCOPY(ngrid,parea,1,area,1)
	725	totarea=SSUM(ngridmx,area,1)
	726
	727	! in 'comdiurn.h'
	728	DO ig=1,ngrid
	729	sinlat(ig)=sin(plat(ig))
	730	coslat(ig)=cos(plat(ig))
	731	sinlon(ig)=sin(plon(ig))
	732	coslon(ig)=cos(plon(ig))
	733	ENDDO
	734
	735	pi=2.*asin(1.) ! NB: pi is a common in comcstfi.h
	736
	737	! managing the tracers, and tests:
	738	! -------------------------------
	739	! Ehouarn: removed; as these tests are now done in initracer.F
	740	! if(tracer) then
	741	!
	742	!! when photochem is used, nqchem_min is the rank
	743	!! of the first chemical species
	744	!
	745	!! Ehouarn: nqchem_min is now meaningless and no longer used
	746	!! nqchem_min = 1
	747	! if (photochem .or. callthermos) then
	748	! chem = .true.
	749	! end if
	750	!
	751	! if (water .or. thermoswater) h2o = .true.
	752	!
	753	!! TESTS
	754	!
	755	! print*,'inifis: TRACERS:'
	756	! write(,) " chem=",chem," h2o=",h2o
	757	!! write(,) " doubleq=",doubleq
	758	!! write(,) " dustbin=",dustbin
	759	!
	760	! if ((doubleq).and.(h2o).and.
	761	! $ (chem)) then
	762	! print*,' 2 dust tracers (doubleq)'
	763	! print*,' 1 water vapour tracer'
	764	! print*,' 1 water ice tracer'
	765	! print*,nqmx-4,' chemistry tracers'
	766	! endif
	767	!
	768	! if ((doubleq).and.(h2o).and.
	769	! $ .not.(chem)) then
	770	! print*,' 2 dust tracers (doubleq)'
	771	! print*,' 1 water vapour tracer'
	772	! print*,' 1 water ice tracer'
	773	! if (nqmx.LT.4) then
	774	! print*,'nqmx should be at least equal to'
	775	! print*,'4 with these options.'
	776	! stop
	777	! endif
	778	! endif
	779	!
	780	! if (.not.(doubleq).and.(h2o).and.
	781	! $ (chem)) then
	782	! if (dustbin.gt.0) then
	783	! print*,dustbin,' dust bins'
	784	! endif
	785	! print*,nqmx-2-dustbin,' chemistry tracers'
	786	! print*,' 1 water vapour tracer'
	787	! print*,' 1 water ice tracer'
	788	! endif
	789	!
	790	! if (.not.(doubleq).and.(h2o).and.
	791	! $ .not.(chem)) then
	792	! if (dustbin.gt.0) then
	793	! print*,dustbin,' dust bins'
	794	! endif
	795	! print*,' 1 water vapour tracer'
	796	! print*,' 1 water ice tracer'
	797	! if (nqmx.gt.(dustbin+2)) then
	798	! print*,'nqmx should be ',(dustbin+2),
	799	! $ ' with these options...'
	800	! print*,'(or check callphys.def)'
	801	! endif
	802	! if (nqmx.lt.(dustbin+2)) then
	803	! write(,) "inifis: nqmx.lt.(dustbin+2)"
	804	! stop
	805	! endif
	806	! endif
	807	!
	808	! endif ! of if (tracer)
	809	!
	810	! RETURN
	811	END

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

Download in other formats:

Original Format