Context Navigation

inifis.F @ 705

Last change on this file since 705 was 705, checked in by emillour, 12 years ago

Mars GCM:

Added possibility to run with a varying EUV cycle following real one. The flag solvarmod=1 triggers this behaviour, with companion flag solvaryear=## , where ## is the Mars Year (from 23 to 30). Setting solvarmod=0 reverts to 'old' behaviour, where there is a constant EUV forcing throughout the run, set by the "solarcondate" flag.
Needs corresponding input data files ("param_v6" subdirectory of "EUV" subdirectory in "datadir").
Added files jthermcalc_e107.F and param_read_e107.F in "aeronomars", modified files euvheat.F90, hrtherm.F, chemthermos.F90, param_v4.h and param_read.F in "aeronomars" and files inifis.F, physiq.F and callkeys.h in "phymars".

FGG

File size: 27.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,",
[677]	181	& "=6 cold (low dust) scenario; =7 warm (high dust) scenario ",
	182	& "=24,25 ... 30 :Mars Year 24, ... or 30 from TES assimilation"
[42]	183	iaervar=3 ! default value
	184	call getin("iaervar",iaervar)
	185	write(,) " iaervar = ",iaervar
	186
[627]	187	write(,) "Reference (visible) dust opacity at 610 Pa ",
[42]	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)"
[577]	251	#ifdef MESOSCALE
	252	nircorr=0 !default value. this is OK below 60 km.
	253	#else
[633]	254	nircorr=0 !default value
[577]	255	#endif
[414]	256	call getin("nircorr",nircorr)
	257	write(,) " nircorr = ",nircorr
	258
[42]	259	write(,) "call turbulent vertical diffusion ?"
	260	calldifv=.true. ! default value
	261	call getin("calldifv",calldifv)
	262	write(,) " calldifv = ",calldifv
	263
[162]	264	write(,) "call thermals ?"
	265	calltherm=.false. ! default value
	266	call getin("calltherm",calltherm)
	267	write(,) " calltherm = ",calltherm
	268
	269	write(,) "output thermal diagnostics ?"
	270	outptherm=.false. ! default value
	271	call getin("outptherm",outptherm)
	272	write(,) " outptherm = ",outptherm
	273
[42]	274	write(,) "call convective adjustment ?"
	275	calladj=.true. ! default value
	276	call getin("calladj",calladj)
	277	write(,) " calladj = ",calladj
	278
[162]	279	if (calltherm .and. (.not. calladj)) then
	280	print*,'Convadj has to be activated when using thermals'
	281	stop
	282	endif
[42]	283
[284]	284	write(,) "call Richardson-based surface layer ?"
	285	callrichsl=.false. ! default value
	286	call getin("callrichsl",callrichsl)
	287	write(,) " callrichsl = ",callrichsl
	288
[288]	289	if (calltherm .and. .not.callrichsl) then
	290	print*,'WARNING WARNING WARNING'
	291	print*,'if calltherm=T we strongly advise that '
	292	print*,'you use the new surface layer scheme '
	293	print*,'by setting callrichsl=T '
	294	endif
	295
[553]	296	if (calladj .and. callrichsl .and. (.not. calltherm)) then
[551]	297	print*,'You should not be calling the convective adjustment
[553]	298	& scheme with the Richardson surface-layer and without the thermals
	299	&. This approach is not
[551]	300	& physically consistent and can lead to unrealistic friction
	301	& values.'
	302	print*,'If you want to use the Ri. surface-layer, either
	303	& activate thermals OR de-activate the convective adjustment.'
	304	stop
	305	endif
[566]	306
[42]	307	write(,) "call CO2 condensation ?"
	308	callcond=.true. ! default value
	309	call getin("callcond",callcond)
	310	write(,) " callcond = ",callcond
	311
	312	write(,)"call thermal conduction in the soil ?"
	313	callsoil=.true. ! default value
	314	call getin("callsoil",callsoil)
	315	write(,) " callsoil = ",callsoil
	316
	317
	318	write(,)"call Lott's gravity wave/subgrid topography ",
	319	& "scheme ?"
	320	calllott=.true. ! default value
	321	call getin("calllott",calllott)
	322	write(,)" calllott = ",calllott
	323
	324
	325	write(,)"rad.transfer is computed every iradia",
	326	& " physical timestep"
	327	iradia=1 ! default value
	328	call getin("iradia",iradia)
	329	write(,)" iradia = ",iradia
	330
	331
	332	write(,)"Output of the exchange coefficient mattrix ?",
	333	& "(for diagnostics only)"
	334	callg2d=.false. ! default value
	335	call getin("callg2d",callg2d)
	336	write(,)" callg2d = ",callg2d
	337
	338	write(,)"Rayleigh scattering : (should be .false. for now)"
	339	rayleigh=.false.
	340	call getin("rayleigh",rayleigh)
	341	write(,)" rayleigh = ",rayleigh
	342
	343
	344	! TRACERS:
	345
[358]	346	! dustbin
[42]	347	write(,)"Transported dust ? (if >0, use 'dustbin' dust bins)"
	348	dustbin=0 ! default value
	349	call getin("dustbin",dustbin)
	350	write(,)" dustbin = ",dustbin
[358]	351	! active
[42]	352	write(,)"Radiatively active dust ? (matters if dustbin>0)"
	353	active=.false. ! default value
	354	call getin("active",active)
	355	write(,)" active = ",active
	356
	357	! Test of incompatibility:
	358	! if active is used, then dustbin should be > 0
	359
	360	if (active.and.(dustbin.lt.1)) then
	361	print*,'if active is used, then dustbin should > 0'
	362	stop
	363	endif
[358]	364	! doubleq
[42]	365	write(,)"use mass and number mixing ratios to predict",
	366	& " dust size ?"
	367	doubleq=.false. ! default value
	368	call getin("doubleq",doubleq)
	369	write(,)" doubleq = ",doubleq
[358]	370	! submicron
[42]	371	submicron=.false. ! default value
	372	call getin("submicron",submicron)
	373	write(,)" submicron = ",submicron
	374
	375	! Test of incompatibility:
	376	! if doubleq is used, then dustbin should be 2
	377
	378	if (doubleq.and.(dustbin.ne.2)) then
	379	print*,'if doubleq is used, then dustbin should be 2'
	380	stop
	381	endif
	382	if (doubleq.and.submicron.and.(nqmx.LT.3)) then
	383	print*,'If doubleq is used with a submicron tracer,'
	384	print*,' then the number of tracers has to be'
	385	print*,' larger than 3.'
	386	stop
	387	endif
[358]	388	! lifting
[42]	389	write(,)"dust lifted by GCM surface winds ?"
	390	lifting=.false. ! default value
	391	call getin("lifting",lifting)
	392	write(,)" lifting = ",lifting
	393
	394	! Test of incompatibility:
	395	! if lifting is used, then dustbin should be > 0
	396
	397	if (lifting.and.(dustbin.lt.1)) then
	398	print*,'if lifting is used, then dustbin should > 0'
	399	stop
	400	endif
[358]	401	! callddevil
[42]	402	write(,)" dust lifted by dust devils ?"
	403	callddevil=.false. !default value
	404	call getin("callddevil",callddevil)
	405	write(,)" callddevil = ",callddevil
	406
	407	! Test of incompatibility:
	408	! if dustdevil is used, then dustbin should be > 0
	409
	410	if (callddevil.and.(dustbin.lt.1)) then
	411	print*,'if dustdevil is used, then dustbin should > 0'
	412	stop
	413	endif
[358]	414	! sedimentation
[42]	415	write(,) "Gravitationnal sedimentation ?"
	416	sedimentation=.true. ! default value
	417	call getin("sedimentation",sedimentation)
	418	write(,) " sedimentation = ",sedimentation
[358]	419	! activice
[42]	420	write(,) "Radiatively active transported atmospheric ",
	421	& "water ice ?"
	422	activice=.false. ! default value
	423	call getin("activice",activice)
	424	write(,) " activice = ",activice
[358]	425	! water
[42]	426	write(,) "Compute water cycle ?"
	427	water=.false. ! default value
	428	call getin("water",water)
	429	write(,) " water = ",water
	430
	431	! Test of incompatibility:
	432
	433	if (activice.and..not.water) then
	434	print*,'if activice is used, water should be used too'
	435	stop
	436	endif
	437
	438	if (water.and..not.tracer) then
	439	print*,'if water is used, tracer should be used too'
	440	stop
	441	endif
[420]	442
[455]	443	! water ice clouds effective variance distribution for sedimentaion
	444	write(,) "effective variance for water ice clouds ?"
	445	nuice_sed=0.45
	446	call getin("nuice_sed",nuice_sed)
	447	write(,) "water_param nueff Sedimentation:", nuice_sed
	448
[420]	449	! ccn factor if no scavenging
	450	write(,) "water param CCN reduc. factor ?", ccn_factor
	451	ccn_factor = 4.5
	452	call getin("ccn_factor",ccn_factor)
	453	write(,)" ccn_factor = ",ccn_factor
	454	write(,)"Careful: only used when microphys=F, otherwise"
	455	write(,)"the contact parameter is used instead;"
	456
[358]	457	! microphys
	458	write(,)"Microphysical scheme for water-ice clouds?"
	459	microphys=.false. ! default value
	460	call getin("microphys",microphys)
	461	write(,)" microphys = ",microphys
[42]	462
[358]	463	! microphysical parameter contact
	464	write(,) "water contact parameter ?"
	465	mteta = 0.95
	466	call getin("mteta",mteta)
	467	write(,) "mteta = ", mteta
	468
	469	! scavenging
	470	write(,)"Dust scavenging by H2O/CO2 snowfall ?"
	471	scavenging=.false. ! default value
	472	call getin("scavenging",scavenging)
	473	write(,)" scavenging = ",scavenging
	474
	475
[42]	476	! Test of incompatibility:
[358]	477	! if scavenging is used, then dustbin should be > 0
[42]	478
[358]	479	if (scavenging.and.(dustbin.lt.1)) then
	480	print*,'if scavenging is used, then dustbin should > 0'
	481	stop
	482	endif
	483	if ((microphys.and..not.doubleq).or.
	484	& (microphys.and..not.active).or.
	485	& (microphys.and..not.scavenging).or.
	486	& (microphys.and..not.water)) then
	487	print*,'if microphys is used, then doubleq,'
	488	print*,'active, water and scavenging must all be used!'
	489	stop
	490	endif
	491	if ((scavenging.and..not.doubleq).or.
	492	& (scavenging.and..not.active).or.
	493	& (scavenging.and..not.water).or.
	494	& (scavenging.and..not.microphys)) then
	495	print*,'if scavenging is used, then doubleq,'
	496	print*,'active, water and microphys must all be used!'
	497	stop
	498	endif
	499
	500	! Test of incompatibility:
	501
[42]	502	write(,) "Permanent water caps at poles ?",
	503	& " .true. is RECOMMENDED"
	504	write(,) "(with .true., North cap is a source of water ",
	505	& "and South pole is a cold trap)"
	506	caps=.true. ! default value
	507	call getin("caps",caps)
	508	write(,) " caps = ",caps
	509
[358]	510	! albedo_h2o_ice
[283]	511	write(,) "water ice albedo ?"
	512	albedo_h2o_ice=0.45
	513	call getin("albedo_h2o_ice",albedo_h2o_ice)
	514	write(,) " albedo_h2o_ice = ",albedo_h2o_ice
[358]	515	! inert_h2o_ice
[283]	516	write(,) "water ice thermal inertia ?"
	517	inert_h2o_ice=2400 ! (J.m^-2.K^-1.s^-1/2)
	518	call getin("inert_h2o_ice",inert_h2o_ice)
	519	write(,) " inert_h2o_ice = ",inert_h2o_ice
[358]	520	! frost_albedo_threshold
[283]	521	write(,) "frost thickness threshold for albedo ?"
[285]	522	frost_albedo_threshold=0.005 ! 5.4 mic (i.e 0.005 kg.m-2)
[283]	523	call getin("frost_albedo_threshold",
	524	& frost_albedo_threshold)
	525	write(,) " frost_albedo_threshold = ",
	526	& frost_albedo_threshold
	527
[485]	528	! call Titus crocus line -- DEFAULT IS NONE
	529	write(,) "Titus crocus line ?"
	530	tituscap=.false. ! default value
	531	call getin("tituscap",tituscap)
	532	write(,) "tituscap",tituscap
[633]	533
[485]	534
[42]	535	write(,) "photochemistry: include chemical species"
	536	photochem=.false. ! default value
	537	call getin("photochem",photochem)
	538	write(,) " photochem = ",photochem
	539
	540
	541	! THERMOSPHERE
	542
	543	write(,) "call thermosphere ?"
	544	callthermos=.false. ! default value
	545	call getin("callthermos",callthermos)
	546	write(,) " callthermos = ",callthermos
	547
	548
	549	write(,) " water included without cycle ",
	550	& "(only if water=.false.)"
	551	thermoswater=.false. ! default value
	552	call getin("thermoswater",thermoswater)
	553	write(,) " thermoswater = ",thermoswater
	554
	555	write(,) "call thermal conduction ?",
	556	& " (only if callthermos=.true.)"
	557	callconduct=.false. ! default value
	558	call getin("callconduct",callconduct)
	559	write(,) " callconduct = ",callconduct
	560
	561	write(,) "call EUV heating ?",
	562	& " (only if callthermos=.true.)"
	563	calleuv=.false. ! default value
	564	call getin("calleuv",calleuv)
	565	write(,) " calleuv = ",calleuv
	566
	567	write(,) "call molecular viscosity ?",
	568	& " (only if callthermos=.true.)"
	569	callmolvis=.false. ! default value
	570	call getin("callmolvis",callmolvis)
	571	write(,) " callmolvis = ",callmolvis
	572
	573	write(,) "call molecular diffusion ?",
	574	& " (only if callthermos=.true.)"
	575	callmoldiff=.false. ! default value
	576	call getin("callmoldiff",callmoldiff)
	577	write(,) " callmoldiff = ",callmoldiff
	578
	579
	580	write(,) "call thermospheric photochemistry ?",
	581	& " (only if callthermos=.true.)"
	582	thermochem=.false. ! default value
	583	call getin("thermochem",thermochem)
	584	write(,) " thermochem = ",thermochem
	585
[705]	586	write(,) "Method to include solar variability"
	587	write(,) "0-> old method (using solarcondate); ",
	588	& "1-> variability wit E10.7"
	589	solvarmod=1
	590	call getin("solvarmod",solvarmod)
	591	write(,) " solvarmod = ",solvarmod
	592
[42]	593	write(,) "date for solar flux calculation:",
[705]	594	& " (1985 < date < 2002)",
	595	$ " (Only used if solvarmod=0)"
[42]	596	write(,) "(Solar min=1996.4 ave=1993.4 max=1990.6)"
	597	solarcondate=1993.4 ! default value
	598	call getin("solarcondate",solarcondate)
	599	write(,) " solarcondate = ",solarcondate
	600
[705]	601	write(,) "Solar variability as observed for MY: "
	602	write(,) "Only if solvarmod=1"
	603	solvaryear=24
	604	call getin("solvaryear",solvaryear)
	605	write(,) " solvaryear = ",solvaryear
	606
[552]	607	write(,) "UV heating efficiency:",
	608	& "measured values between 0.19 and 0.23 (Fox et al. 1996)",
	609	& "lower values may be used to compensate low 15 um cooling"
	610	euveff=0.21 !default value
	611	call getin("euveff",euveff)
	612	write(,) " euveff = ", euveff
[42]	613
	614	if (.not.callthermos) then
	615	if (thermoswater) then
	616	print*,'if thermoswater is set, callthermos must be true'
	617	stop
	618	endif
	619	if (callconduct) then
	620	print*,'if callconduct is set, callthermos must be true'
	621	stop
	622	endif
	623	if (calleuv) then
	624	print*,'if calleuv is set, callthermos must be true'
	625	stop
	626	endif
	627	if (callmolvis) then
	628	print*,'if callmolvis is set, callthermos must be true'
	629	stop
	630	endif
	631	if (callmoldiff) then
	632	print*,'if callmoldiff is set, callthermos must be true'
	633	stop
	634	endif
	635	if (thermochem) then
	636	print*,'if thermochem is set, callthermos must be true'
	637	stop
	638	endif
	639	endif
	640
	641	! Test of incompatibility:
	642	! if photochem is used, then water should be used too
	643
	644	if (photochem.and..not.water) then
	645	print*,'if photochem is used, water should be used too'
	646	stop
	647	endif
	648
	649	! if callthermos is used, then thermoswater should be used too
	650	! (if water not used already)
	651
	652	if (callthermos .and. .not.water) then
	653	if (callthermos .and. .not.thermoswater) then
	654	print*,'if callthermos is used, water or thermoswater
	655	& should be used too'
	656	stop
	657	endif
	658	endif
	659
	660	PRINT*,'--------------------------------------------'
	661	PRINT*
	662	PRINT*
	663	ELSE
	664	write(,)
	665	write(,) 'Cannot read file callphys.def. Is it here ?'
	666	stop
	667	ENDIF
	668
	669	8000 FORMAT(t5,a12,l8)
	670	8001 FORMAT(t5,a12,i8)
	671
	672	PRINT*
	673	PRINT*,'inifis: daysec',daysec
	674	PRINT*
	675	PRINT*,'inifis: The radiative transfer is computed:'
	676	PRINT*,' each ',iradia,' physical time-step'
	677	PRINT,' or each ',iradiadtphys,' seconds'
	678	PRINT*
	679	! --------------------------------------------------------------
	680	! Managing the Longwave radiative transfer
	681	! --------------------------------------------------------------
	682
	683	! In most cases, the run just use the following values :
	684	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	685	callemis=.true.
	686	! ilwd=10*int(daysec/dtphys) ! bug before 22/10/01
	687	ilwd=1
	688	ilwn=1 !2
	689	ilwb=1 !2
	690	linear=.true.
	691	ncouche=3
	692	alphan=0.4
	693	semi=0
	694
	695	! BUT people working hard on the LW may want to read them in 'radia.def'
	696	! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	697
	698	OPEN(99,file='radia.def',status='old',form='formatted'
	699	. ,iostat=ierr)
	700	IF(ierr.EQ.0) THEN
	701	write(,) 'inifis: Reading radia.def !!!'
	702	READ(99,fmt='(a)') ch1
	703	READ(99,*) callemis
	704	WRITE(*,8000) ch1,callemis
	705
	706	READ(99,fmt='(a)') ch1
	707	READ(99,*) iradia
	708	WRITE(*,8001) ch1,iradia
	709
	710	READ(99,fmt='(a)') ch1
	711	READ(99,*) ilwd
	712	WRITE(*,8001) ch1,ilwd
	713
	714	READ(99,fmt='(a)') ch1
	715	READ(99,*) ilwn
	716	WRITE(*,8001) ch1,ilwn
	717
	718	READ(99,fmt='(a)') ch1
	719	READ(99,*) linear
	720	WRITE(*,8000) ch1,linear
	721
	722	READ(99,fmt='(a)') ch1
	723	READ(99,*) ncouche
	724	WRITE(*,8001) ch1,ncouche
	725
	726	READ(99,fmt='(a)') ch1
	727	READ(99,*) alphan
	728	WRITE(,) ch1,alphan
	729
	730	READ(99,fmt='(a)') ch1
	731	READ(99,*) ilwb
	732	WRITE(*,8001) ch1,ilwb
	733
	734
	735	READ(99,fmt='(a)') ch1
	736	READ(99,'(l1)') callg2d
	737	WRITE(*,8000) ch1,callg2d
	738
	739	READ(99,fmt='(a)') ch1
	740	READ(99,*) semi
	741	WRITE(,) ch1,semi
	742	end if
	743	CLOSE(99)
	744
	745	!-----------------------------------------------------------------------
	746	! Some more initialization:
	747	! ------------------------
	748
[226]	749	! in 'comgeomfi.h'
[42]	750	CALL SCOPY(ngrid,plon,1,long,1)
	751	CALL SCOPY(ngrid,plat,1,lati,1)
	752	CALL SCOPY(ngrid,parea,1,area,1)
	753	totarea=SSUM(ngridmx,area,1)
	754
	755	! in 'comdiurn.h'
	756	DO ig=1,ngrid
	757	sinlat(ig)=sin(plat(ig))
	758	coslat(ig)=cos(plat(ig))
	759	sinlon(ig)=sin(plon(ig))
	760	coslon(ig)=cos(plon(ig))
	761	ENDDO
	762
	763	pi=2.*asin(1.) ! NB: pi is a common in comcstfi.h
	764
	765	! managing the tracers, and tests:
	766	! -------------------------------
	767	! Ehouarn: removed; as these tests are now done in initracer.F
	768	! if(tracer) then
	769	!
	770	!! when photochem is used, nqchem_min is the rank
	771	!! of the first chemical species
	772	!
	773	!! Ehouarn: nqchem_min is now meaningless and no longer used
	774	!! nqchem_min = 1
	775	! if (photochem .or. callthermos) then
	776	! chem = .true.
	777	! end if
	778	!
	779	! if (water .or. thermoswater) h2o = .true.
	780	!
	781	!! TESTS
	782	!
	783	! print*,'inifis: TRACERS:'
	784	! write(,) " chem=",chem," h2o=",h2o
	785	!! write(,) " doubleq=",doubleq
	786	!! write(,) " dustbin=",dustbin
	787	!
	788	! if ((doubleq).and.(h2o).and.
	789	! $ (chem)) then
	790	! print*,' 2 dust tracers (doubleq)'
	791	! print*,' 1 water vapour tracer'
	792	! print*,' 1 water ice tracer'
	793	! print*,nqmx-4,' chemistry tracers'
	794	! endif
	795	!
	796	! if ((doubleq).and.(h2o).and.
	797	! $ .not.(chem)) then
	798	! print*,' 2 dust tracers (doubleq)'
	799	! print*,' 1 water vapour tracer'
	800	! print*,' 1 water ice tracer'
	801	! if (nqmx.LT.4) then
	802	! print*,'nqmx should be at least equal to'
	803	! print*,'4 with these options.'
	804	! stop
	805	! endif
	806	! endif
	807	!
	808	! if (.not.(doubleq).and.(h2o).and.
	809	! $ (chem)) then
	810	! if (dustbin.gt.0) then
	811	! print*,dustbin,' dust bins'
	812	! endif
	813	! print*,nqmx-2-dustbin,' chemistry tracers'
	814	! print*,' 1 water vapour tracer'
	815	! print*,' 1 water ice tracer'
	816	! endif
	817	!
	818	! if (.not.(doubleq).and.(h2o).and.
	819	! $ .not.(chem)) then
	820	! if (dustbin.gt.0) then
	821	! print*,dustbin,' dust bins'
	822	! endif
	823	! print*,' 1 water vapour tracer'
	824	! print*,' 1 water ice tracer'
	825	! if (nqmx.gt.(dustbin+2)) then
	826	! print*,'nqmx should be ',(dustbin+2),
	827	! $ ' with these options...'
	828	! print*,'(or check callphys.def)'
	829	! endif
	830	! if (nqmx.lt.(dustbin+2)) then
	831	! write(,) "inifis: nqmx.lt.(dustbin+2)"
	832	! stop
	833	! endif
	834	! endif
	835	!
	836	! endif ! of if (tracer)
	837	!
	838	! RETURN
	839	END

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

Download in other formats:

Original Format