Context Navigation

rcm1d.F @ 3982

Last change on this file since 3982 was 3978, checked in by debatzbr, 3 months ago
Pluto PCM: Initialize gas profiles for clouds in 1D. BBT
Property svn:executable set to ``*
File size: 39.0 KB

Rev	Line
[3184]	1	program rcm1d
	2
	3	! to use 'getin'
	4	use ioipsl_getincom, only: getin
	5	use dimphy, only : init_dimphy
	6	use mod_grid_phy_lmdz, only : regular_lonlat
	7	use infotrac, only: nqtot, tname
[3585]	8	use tracer_h, only: noms
[3184]	9	use surfdat_h, only: albedodat, phisfi, dryness,
	10	& zmea, zstd, zsig, zgam, zthe,
	11	& emissiv, emisice, iceradius,
[3916]	12	& dtemisice, n2frac
[3184]	13	use comdiurn_h, only: sinlat, coslat, sinlon, coslon
	14	use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat, volcapa
	15	use phyredem, only: physdem0,physdem1
[3816]	16	use geometry_mod, only: init_geometry
[3184]	17	use planete_mod, only: apoastr,periastr,year_day,peri_day,
[3542]	18	& obliquit,z0,lmixmin,emin_turb,coefvis,coefir,
[3184]	19	& timeperi,e_elips,p_elips
	20	use comcstfi_mod, only: pi, cpp, rad, g, r,
	21	& mugaz, rcp, omeg
[3749]	22	use time_phylmdz_mod, only: daysec, dtphys, diagfi_output_rate,
[3758]	23	& nday
[3978]	24	use callkeys_mod, only: tracer, specOLR,pceil,haze,
	25	& callmufi,callmuclouds
[3184]	26	USE comvert_mod, ONLY: ap,bp,aps,bps,pa,preff, sig,
	27	& presnivs,pseudoalt,scaleheight
	28	USE vertical_layers_mod, ONLY: init_vertical_layers
	29	USE logic_mod, ONLY: hybrid
	30	use radinc_h, only: naerkind
	31	use regular_lonlat_mod, only: init_regular_lonlat
	32	use planete_mod, only: ini_planete_mod
	33	use physics_distribution_mod, only: init_physics_distribution
	34	use mod_interface_dyn_phys, only: init_interface_dyn_phys
	35	use inifis_mod, only: inifis
[3672]	36	use phys_state_var_mod, only: phys_state_var_init, tsurf, tsoil
[3184]	37	use physiq_mod, only: physiq
	38	implicit none
	39
	40	!==================================================================
[3258]	41	!
[3184]	42	! Purpose
	43	! -------
	44	! Run the physics package of the universal model in a 1D column.
[3258]	45	!
[3184]	46	! It can be compiled with a command like (e.g. for 25 layers):
	47	! "makegcm -p std -d 25 rcm1d"
	48	! It requires the files "callphys.def", "z2sig.def",
	49	! "traceur.def", and "run.def" with a line "INCLUDEDEF=callphys.def"
	50	!
	51	! Authors
	52	! -------
	53	! Frederic Hourdin
	54	! R. Fournier
	55	! F. Forget
	56	! F. Montmessin (water ice added)
	57	! R. Wordsworth
	58	! B. Charnay
	59	! A. Spiga
	60	! J. Leconte (2012)
	61	!
	62	!==================================================================
	63
	64	#include "dimensions.h"
	65	#include "paramet.h"
	66	!include "dimphys.h"
	67	#include "netcdf.inc"
	68	#include "comgeom.h"
	69
	70	c --------------------------------------------------------------
	71	c Declarations
	72	c --------------------------------------------------------------
	73	c
	74	INTEGER unitstart ! unite d'ecriture de "startfi"
	75	INTEGER nlayer,nlevel,nsoil,ndt
	76	INTEGER ilayer,ilevel,isoil,idt,iq
	77	LOGICAl firstcall,lastcall
	78	c
	79	INTEGER day0 ! date initial (sol ; =0 a Ls=0)
[3329]	80	INTEGER lecttsoil ! lecture of tsoil from proftsoil
	81	INTEGER lecthaze ! lecture of haze from profhaze
[3978]	82	INTEGER lectC2H2 ! lecture of gases from profC2H2
	83	INTEGER lectC2H6 ! lecture of gases from profC2H6
	84	INTEGER lectC4H2 ! lecture of gases from profC4H2
	85	INTEGER lectC6H6 ! lecture of gases from profC6H6
	86	INTEGER lectHCN ! lecture of gases from profHCN
[3749]	87	REAL day ! date during the run
	88	INTEGER day_step ! number of time steps per day
[3184]	89	REAL time ! time (0<time<1 ; time=0.5 a midi)
	90	REAL play(llm) ! Pressure at the middle of the layers (Pa)
	91	REAL plev(llm+1) ! intermediate pressure levels (pa)
[3672]	92	REAL psurf
[3184]	93	REAL u(llm),v(llm) ! zonal, meridional wind
	94	REAL gru,grv ! prescribed "geostrophic" background wind
	95	REAL temp(llm) ! temperature at the middle of the layers
	96	REAL,ALLOCATABLE :: q(:,:) ! tracer mixing ratio (e.g. kg/kg)
	97	REAL,ALLOCATABLE :: qsurf(:) ! tracer surface budget (e.g. kg.m-2)
[3232]	98	! REAL n2ice ! n2ice layer (kg.m-2) !not used anymore
[3258]	99	integer :: i_n2=0 ! tracer index of n2 ice
	100	integer :: i_ch4_ice=0 ! tracer index of ch4 ice
	101	integer :: i_ch4_gas=0 ! tracer index of ch4 gas
	102	integer :: i_co_ice=0 ! tracer index of co ice
	103	integer :: i_co_gas=0 ! tracer index of co gas
[3978]	104	integer :: i_C2H2_mugas=0 ! tracer index of C2H2 gas
	105	integer :: i_C2H6_mugas=0 ! tracer index of C2H6 gas
	106	integer :: i_C4H2_mugas=0 ! tracer index of C4H2 gas
	107	integer :: i_C6H6_mugas=0 ! tracer index of C6H6 gas
	108	integer :: i_HCN_mugas=0 ! tracer index of HCN gas
[3258]	109	integer :: i_prec_haze=0 ! tracer index of haze
	110	integer :: i_haze=0 ! tracer index of haze
	111	integer :: i_haze10=0 ! tracer index of haze
	112	integer :: i_haze30=0 ! tracer index of haze
	113	integer :: i_haze50=0 ! tracer index of haze
	114	integer :: i_haze100=0 ! tracer index of haze
	115
[3184]	116	REAL emis(1) ! surface layer
	117	REAL q2(llm+1) ! Turbulent Kinetic Energy
	118	REAL zlay(llm) ! altitude estimee dans les couches (km)
	119
	120	c Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s)
	121	REAL du(llm),dv(llm),dtemp(llm)
	122	REAL dudyn(llm),dvdyn(llm),dtempdyn(llm)
[3258]	123	REAL dpsurf(1)
[3184]	124	REAL,ALLOCATABLE :: dq(:,:)
	125	REAL,ALLOCATABLE :: dqdyn(:,:)
	126
	127	c Various intermediate variables
	128	! INTEGER thermo
	129	REAL zls
	130	REAL phi(llm),h(llm),s(llm)
	131	REAL pks, ptif, w(llm)
	132	INTEGER ierr, aslun
	133	REAL tmp1(0:llm),tmp2(0:llm)
	134	integer :: nq !=1 ! number of tracers
[3258]	135
[3184]	136	character*2 str2
	137	character (len=7) :: str7
	138	character(len=44) :: txt
	139	character(len=44) :: tracer_profile_file_name
	140
	141	logical oldcompare, earthhack,saveprofile
	142
	143	! added by RW for zlay computation
	144	real Hscale, Hmax, rho, dz
	145
[3258]	146	! pluto specific
	147	real ch4ref ! % ch4
	148	real coref ! % co
	149	real hazeref ! haze kg/kg
	150	real prechazeref ! prec haze kg/kg
	151
	152
[3184]	153	! added by RW for autozlevs computation
	154	logical autozlevs
	155	real nu, xx, pMIN, zlev, Htop
	156	real logplevs(llm)
	157
	158	! added by BC to accelerate convergence
	159	logical accelerate_temperature
	160	real factor_acceleration
	161
	162	! added by AS to avoid the use of adv trac common
	163	character*30,allocatable :: nametmp(:) !
	164
	165	real :: latitude(1), longitude(1), cell_area(1)
	166
	167	! added by JVO and YJ to read modern traceur.def
	168	character(len=500) :: line ! to store a line of text
	169	LOGICAL :: moderntracdef=.false. ! JVO, YJ : modern traceur.def
	170
	171	c=======================================================================
	172	c INITIALISATION
	173	c=======================================================================
[3258]	174	! check if 'rcm1d.def' file is around
[3184]	175	open(90,file='rcm1d.def',status='old',form='formatted',
	176	& iostat=ierr)
	177	if (ierr.ne.0) then
	178	write(,) 'Cannot find required file "rcm1d.def"'
	179	write(,) 'which should contain some input parameters'
	180	write(,) ' ... might as well stop here ...'
	181	stop
	182	else
	183	close(90)
	184	endif
	185
	186	! now, run.def is needed anyway. so we create a dummy temporary one
	187	! for ioipsl to work. if a run.def is already here, stop the
	188	! program and ask the user to do a bit of cleaning
	189	open(90,file='run.def',status='old',form='formatted',
	190	& iostat=ierr)
	191	if (ierr.eq.0) then
	192	close(90)
	193	write(,) 'There is already a run.def file.'
	194	write(,) 'This is not compatible with 1D runs.'
	195	write(,) 'Please remove the file and restart the run.'
	196	write(,) 'Runtime parameters are supposed to be in rcm1d.def'
	197	stop
	198	else
	199	call system('touch run.def')
	200	call system("echo 'INCLUDEDEF=callphys.def' >> run.def")
	201	call system("echo 'INCLUDEDEF=rcm1d.def' >> run.def")
	202	endif
	203
	204	! read nq from traceur.def
	205	open(90,file='traceur.def',status='old',form='formatted',
	206	& iostat=ierr)
	207	if (ierr.eq.0) then
	208	! read number of tracers:
	209	!! - Modif. by JVO and YJ for modern planetary traceur.def ---------------
	210	READ(90,'(A)') line
	211	IF (trim(line).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def
	212	READ(line,*,iostat=ierr) nq ! Try standard traceur.def
	213	ELSE
	214	moderntracdef = .true.
	215	DO
	216	READ(90,'(A)',iostat=ierr) line
	217	IF (ierr.eq.0) THEN
	218	IF (index(line,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
	219	READ(line,*,iostat=ierr) nq
	220	EXIT
	221	ENDIF
	222	ELSE ! If pb, or if reached EOF without having found nbtr
	223	write(,) "rcm1d: error reading number of tracers"
	224	write(,) " (first line of traceur.def) "
	225	stop
	226	ENDIF
	227	ENDDO
	228	ENDIF ! if modern or standard traceur.def
	229	else
	230	nq=0
	231	endif
	232	close(90)
[3258]	233
[3184]	234	! Initialize dimphy module
[3258]	235	call init_dimphy(1,llm)
	236
[3184]	237	! now initialize arrays using phys_state_var_init
	238	! but first initialise naerkind (from callphys.def)
	239	naerkind=0 !default
	240	call getin("naerkind",naerkind)
[3258]	241
[3184]	242	call phys_state_var_init(nq)
[3258]	243
[3184]	244	saveprofile=.false.
	245	saveprofile=.true.
	246
	247	c ----------------------------------------
	248	c Default values (corresponding to Mars)
	249	c ----------------------------------------
	250
	251	pi=2.E+0*asin(1.E+0)
	252
[3258]	253	c Parametres Couche limite et Turbulence
[3184]	254	c --------------------------------------
[3258]	255	z0 = 1.e-2 ! surface roughness (m) ~0.01
[3184]	256	emin_turb = 1.e-6 ! energie minimale ~1.e-8
	257	lmixmin = 30 ! longueur de melange ~100
[3258]	258
[3184]	259	c propriete optiques des calottes et emissivite du sol
	260	c ----------------------------------------------------
	261	emissiv= 0.95 ! Emissivite du sol martien ~.95
	262	emisice(1)=0.95 ! Emissivite calotte nord
[3258]	263	emisice(2)=0.95 ! Emissivite calotte sud
[3184]	264
[3232]	265	iceradius(1) = 100.e-6 ! mean scat radius of n2 snow (north)
	266	iceradius(2) = 100.e-6 ! mean scat radius of n2 snow (south)
[3184]	267	dtemisice(1) = 2. ! time scale for snow metamorphism (north)
	268	dtemisice(2) = 2. ! time scale for snow metamorphism (south
	269	hybrid=.false.
	270
	271	c ------------------------------------------------------
[3258]	272	c Load parameters from "run.def" and "gases.def"
[3184]	273	c ------------------------------------------------------
	274
	275
	276	! check if we are going to run with or without tracers
	277	write(,) "Run with or without tracer transport ?"
	278	tracer=.false. ! default value
	279	call getin("tracer",tracer)
	280	write(,) " tracer = ",tracer
	281
	282	! OK. now that run.def has been read once -- any variable is in memory.
	283	! so we can dump the dummy run.def
	284	! call system("rm -rf run.def") ! Ehouarn: delay this to after inifis
	285
	286	! while we're at it, check if there is a 'traceur.def' file
	287	! and preocess it, if necessary. Otherwise initialize tracer names
	288	if (tracer) then
	289	! load tracer names from file 'traceur.def'
	290	open(90,file='traceur.def',status='old',form='formatted',
	291	& iostat=ierr)
	292	if (ierr.eq.0) then
	293	write(,) "rcm1d: Reading file traceur.def"
	294	! read number of tracers:
	295	!! - Modif. by JVO and YJ for modern planetary traceur.def ---------------
	296	READ(90,'(A)') line
	297	IF (trim(line).ne.'#ModernTrac-v1') THEN ! Test modern traceur.def
	298	READ(line,*,iostat=ierr) nq ! Try standard traceur.def
	299	ELSE
	300	moderntracdef = .true.
	301	DO
	302	READ(90,'(A)',iostat=ierr) line
	303	IF (ierr.eq.0) THEN
	304	IF (index(line,'#').ne.1) THEN ! Allows arbitary number of comments lines in the header
	305	READ(line,*,iostat=ierr) nq
	306	EXIT
	307	ENDIF
	308	ELSE ! If pb, or if reached EOF without having found nbtr
	309	write(,) "rcm1d: error reading number of tracers"
	310	write(,) " (first line of traceur.def) "
	311	stop
	312	ENDIF
	313	ENDDO
	314	ENDIF ! if modern or standard traceur.def
	315	nqtot=nq ! set value of nqtot (in infotrac module) as nq
	316	if (ierr.ne.0) then
	317	write(,) "rcm1d: error reading number of tracers"
	318	write(,) " (first line of traceur.def) "
	319	stop
	320	endif
	321	if (nq>0) then
	322	allocate(tname(nq))
	323	allocate(noms(nq))
	324	allocate(q(llm,nq))
	325	allocate(qsurf(nq))
	326	allocate(dq(llm,nq))
	327	allocate(dqdyn(llm,nq))
	328	else
	329	write(,) "rcm1d: Error. You set tracer=.true."
	330	write(,) " but # of tracers in traceur.def is ",nq
	331	stop
	332	endif
[3258]	333
[3184]	334	do iq=1,nq
	335	! minimal version, just read in the tracer names, 1 per line
	336	read(90,*,iostat=ierr) tname(iq)
	337	noms(iq)=tname(iq)
	338	if (ierr.ne.0) then
	339	write(,) 'rcm1d: error reading tracer names...'
	340	stop
	341	endif
	342	enddo !of do iq=1,nq
[3329]	343	! check for n2 / ch4_ice tracers:
[3258]	344	i_n2=0
[3184]	345	do iq=1,nq
[3232]	346	if (tname(iq)=="n2") then
[3258]	347	i_n2=iq
	348	elseif (tname(iq)=="ch4_ice") then
	349	i_ch4_ice=iq
	350	elseif (tname(iq)=="co_ice") then
	351	i_co_ice=iq
	352	elseif (tname(iq)=="ch4_gas") then
	353	i_ch4_gas=iq
	354	elseif (tname(iq)=="co_gas") then
	355	i_co_gas=iq
[3978]	356	elseif (tname(iq)=="C2H2_mugas") then
	357	i_C2H2_mugas=iq
	358	elseif (tname(iq)=="C2H6_mugas") then
	359	i_C2H6_mugas=iq
	360	elseif (tname(iq)=="C4H2_mugas") then
	361	i_C4H2_mugas=iq
	362	elseif (tname(iq)=="C6H6_mugas") then
	363	i_C6H6_mugas=iq
	364	elseif (tname(iq)=="HCN_mugas") then
	365	i_HCN_mugas=iq
[3329]	366	elseif (tname(iq)=="haze") then
	367	i_haze=iq
[3184]	368	endif
	369	enddo
	370	else
	371	write(,) 'Cannot find required file "traceur.def"'
	372	write(,) ' If you want to run with tracers, I need it'
	373	write(,) ' ... might as well stop here ...'
	374	stop
	375	endif
	376	close(90)
	377
	378
	379	else ! of if (tracer)
	380	nqtot=0
	381	nq=0
	382	! still, make allocations for 1 dummy tracer
[3258]	383	allocate(tname(1))
[3184]	384	allocate(qsurf(1))
	385	allocate(q(llm,1))
	386	allocate(dq(llm,1))
[3258]	387
[3184]	388	! Check that tracer boolean is compliant with number of tracers
	389	! -- otherwise there is an error (and more generally we have to be consistent)
	390	if (nq .ge. 1) then
	391	write(,) "------------------------------"
	392	write(,) "rcm1d: You set tracer=.false."
	393	write(,) " But set number of tracers to ",nq
	394	write(,) " > If you want tracers, set tracer=.true."
	395	write(,) "------------------------------"
	396	stop
	397	endif
	398	endif ! of if (tracer)
	399
	400	!!! GEOGRAPHICAL INITIALIZATIONS
	401	!!! AREA. useless in 1D
	402	cell_area(1)=1.E+0 !JL+EM to have access to the area in the diagfi.nc files.
	403	!!! GEOPOTENTIAL. useless in 1D because control by surface pressure
	404	phisfi(1)=0.E+0
[3258]	405	!!! LATITUDE. can be set.
[3184]	406	latitude=0 ! default value for latitude
	407	PRINT *,'latitude (in degrees) ?'
	408	call getin("latitude",latitude)
	409	write(,) " latitude = ",latitude
	410	latitude=latitude*pi/180.E+0
	411	!!! LONGITUDE. useless in 1D.
	412	longitude=0.E+0
	413	longitude=longitude*pi/180.E+0
	414
	415
	416	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	417	!!!! PLANETARY CONSTANTS !!!!
	418	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	419
	420	g = -99999.
	421	PRINT *,'GRAVITY in m s-2 ?'
	422	call getin("g",g)
	423	IF (g.eq.-99999.) THEN
	424	PRINT *,"STOP. I NEED g IN RCM1D.DEF."
	425	STOP
	426	ELSE
	427	PRINT *,"--> g = ",g
	428	ENDIF
	429
	430	rad = -99999.
	431	PRINT *,'PLANETARY RADIUS in m ?'
	432	call getin("rad",rad)
	433	! Planetary radius is needed to compute shadow of the rings
[3232]	434	IF (rad.eq.-99999.) THEN
[3184]	435	PRINT *,"STOP. I NEED rad IN RCM1D.DEF."
	436	STOP
	437	ELSE
	438	PRINT *,"--> rad = ",rad
	439	ENDIF
	440
	441	daysec = -99999.
	442	PRINT *,'LENGTH OF A DAY in s ?'
	443	call getin("daysec",daysec)
	444	IF (daysec.eq.-99999.) THEN
	445	PRINT *,"STOP. I NEED daysec IN RCM1D.DEF."
	446	STOP
	447	ELSE
	448	PRINT *,"--> daysec = ",daysec
	449	ENDIF
	450	omeg=4.*asin(1.)/(daysec)
	451	PRINT *,"OK. FROM THIS I WORKED OUT:"
	452	PRINT *,"--> omeg = ",omeg
	453
	454	year_day = -99999.
	455	PRINT *,'LENGTH OF A YEAR in days ?'
	456	call getin("year_day",year_day)
	457	IF (year_day.eq.-99999.) THEN
	458	PRINT *,"STOP. I NEED year_day IN RCM1D.DEF."
	459	STOP
[3258]	460	ELSE
[3184]	461	PRINT *,"--> year_day = ",year_day
	462	ENDIF
	463
	464	periastr = -99999.
	465	PRINT *,'MIN DIST STAR-PLANET in AU ?'
	466	call getin("periastr",periastr)
	467	IF (periastr.eq.-99999.) THEN
	468	PRINT *,"STOP. I NEED periastr IN RCM1D.DEF."
	469	STOP
	470	ELSE
	471	PRINT *,"--> periastr = ",periastr
	472	ENDIF
	473
	474	apoastr = -99999.
	475	PRINT *,'MAX DIST STAR-PLANET in AU ?'
	476	call getin("apoastr",apoastr)
	477	IF (apoastr.eq.-99999.) THEN
	478	PRINT *,"STOP. I NEED apoastr IN RCM1D.DEF."
	479	STOP
	480	ELSE
	481	PRINT *,"--> apoastr = ",apoastr
	482	ENDIF
	483
	484	peri_day = -99999.
	485	PRINT *,'DATE OF PERIASTRON in days ?'
	486	call getin("peri_day",peri_day)
	487	IF (peri_day.eq.-99999.) THEN
	488	PRINT *,"STOP. I NEED peri_day IN RCM1D.DEF."
	489	STOP
	490	ELSE IF (peri_day.gt.year_day) THEN
	491	PRINT *,"STOP. peri_day.gt.year_day"
	492	STOP
[3258]	493	ELSE
[3184]	494	PRINT *,"--> peri_day = ", peri_day
[3258]	495	ENDIF
[3184]	496
[3258]	497	obliquit = -99999.
[3184]	498	PRINT *,'OBLIQUITY in deg ?'
	499	call getin("obliquit",obliquit)
	500	IF (obliquit.eq.-99999.) THEN
	501	PRINT *,"STOP. I NEED obliquit IN RCM1D.DEF."
	502	STOP
	503	ELSE
	504	PRINT *,"--> obliquit = ",obliquit
[3258]	505	ENDIF
[3184]	506
	507	psurf = -99999.
	508	PRINT *,'SURFACE PRESSURE in Pa ?'
	509	call getin("psurf",psurf)
	510	IF (psurf.eq.-99999.) THEN
	511	PRINT *,"STOP. I NEED psurf IN RCM1D.DEF."
	512	STOP
	513	ELSE
	514	PRINT *,"--> psurf = ",psurf
	515	ENDIF
	516	!! we need reference pressures
	517	pa=psurf/30.
	518	preff=psurf ! these values are not needed in 1D (are you sure JL12)
	519
	520	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	521	!!!! END PLANETARY CONSTANTS !!!!
	522	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	523
	524	c Date et heure locale du debut du run
	525	c ------------------------------------
	526	c Date (en sols depuis le solstice de printemps) du debut du run
	527	day0 = 0 ! default value for day0
	528	write(,) 'Initial date (in martian sols ; =0 at Ls=0)?'
	529	call getin("day0",day0)
	530	day=float(day0)
	531	write(,) " day0 = ",day0
	532	c Heure de demarrage
	533	time=0 ! default value for time
	534	write(,)'Initial local time (in hours, between 0 and 24)?'
	535	call getin("time",time)
	536	write(,)" time = ",time
	537	time=time/24.E+0 ! convert time (hours) to fraction of sol
	538
	539
	540	c Discretisation (Definition de la grille et des pas de temps)
	541	c --------------
	542	c
	543	nlayer=llm
	544	nlevel=nlayer+1
	545
	546	day_step=48 ! default value for day_step
	547	PRINT *,'Number of time steps per sol ?'
	548	call getin("day_step",day_step)
	549	write(,) " day_step = ",day_step
	550
[3749]	551	diagfi_output_rate=24 ! default value for diagfi_output_rate
[3184]	552	PRINT *,'Nunber of steps between writediagfi ?'
[3749]	553	call getin("diagfi_output_rate",diagfi_output_rate)
	554	write(,) " diagfi_output_rate = ",diagfi_output_rate
[3184]	555
	556	ndt=10 ! default value for ndt
	557	PRINT *,'Number of sols to run ?'
	558	call getin("ndt",ndt)
	559	write(,) " ndt = ",ndt
	560	nday=ndt
	561
[3258]	562	ndt=ndt*day_step
	563	dtphys=daysec/day_step
[3184]	564	write(,)"-------------------------------------"
	565	write(,)"-------------------------------------"
	566	write(,)"--> Physical timestep is ",dtphys
	567	write(,)"-------------------------------------"
	568	write(,)"-------------------------------------"
	569
	570	! initializations, as with iniphysiq.F90 for the 3D GCM
	571	call init_physics_distribution(regular_lonlat,4,
	572	& 1,1,1,nlayer,1)
	573	call init_interface_dyn_phys
	574	CALL init_regular_lonlat(1,1,longitude,latitude,
	575	& (/0.,0./),(/0.,0./))
	576	call init_geometry(1,longitude,latitude,
	577	& (/0.,0.,0.,0./),(/0.,0.,0.,0./),
[3718]	578	& cell_area,[1])
[3184]	579	! Ehouarn: init_vertial_layers called later (because disvert not called yet)
	580	! call init_vertical_layers(nlayer,preff,scaleheight,
	581	! & ap,bp,aps,bps,presnivs,pseudoalt)
	582	! call init_dimphy(1,nlayer) ! Initialize dimphy module
	583	call ini_planete_mod(nlayer,preff,ap,bp)
	584
	585	!!! CALL INIFIS
	586	!!! - read callphys.def
	587	!!! - calculate sine and cosine of longitude and latitude
[3258]	588	!!! - calculate mugaz and cp
[3184]	589	!!! NB: some operations are being done dummily in inifis in 1D
	590	!!! - physical constants: nevermind, things are done allright below
	591	!!! - physical frequency: nevermind, in inifis this is a simple print
	592	cpp=-9999. ! dummy init for inifis, will be rewrite later on
	593	r=-9999. ! dummy init for inifis, will be rewrite later on
	594	CALL inifis(1,llm,nq,day0,daysec,nday,dtphys,
	595	. latitude,longitude,cell_area,rad,g,r,cpp)
	596
	597	nsoil=nsoilmx
	598
	599	! At this point, both getin() and getin_p() functions have been used,
	600	! and the run.def file can be removed.
	601	call system("rm -rf run.def")
	602
	603	!!! We check everything is OK.
	604	PRINT *,"CHECK"
	605	PRINT *,"--> mugaz = ",mugaz
	606	PRINT *,"--> cpp = ",cpp
	607	r = 8.314511E+0 * 1000.E+0 / mugaz
	608	rcp = r / cpp
	609
	610	c output spectrum?
	611	write(,)"Output spectral OLR?"
	612	specOLR=.false.
	613	call getin("specOLR",specOLR)
	614	write(,)" specOLR = ",specOLR
	615
	616	c option for temperature evolution?
	617	write(,)"accelerate temperature?"
	618	accelerate_temperature=.false.
	619	call getin("accelerate_temperature",accelerate_temperature)
	620	write(,)" accelerate_temperature = ",accelerate_temperature
	621
	622	write(,)"factor for temperature acceleration"
	623	factor_acceleration=1.0
	624	call getin("factor_acceleration",factor_acceleration)
	625	write(,)" factor_acceleration = ",factor_acceleration
	626
	627
	628
	629	c
	630	c pour le schema d'ondes de gravite
	631	c ---------------------------------
	632	zmea(1)=0.E+0
	633	zstd(1)=0.E+0
	634	zsig(1)=0.E+0
	635	zgam(1)=0.E+0
	636	zthe(1)=0.E+0
	637
	638	c Initialisation des traceurs
	639	c ---------------------------
	640
	641	DO iq=1,nq
	642	DO ilayer=1,nlayer
	643	q(ilayer,iq) = 0.
	644	ENDDO
	645	ENDDO
[3258]	646
[3184]	647	DO iq=1,nq
	648	qsurf(iq) = 0.
	649	ENDDO
[3258]	650
	651	if (tracer) then ! Initialize tracers here.
	652
[3184]	653	write(,) "rcm1d : initializing tracers profiles"
	654
	655	do iq=1,nq
	656
[3258]	657	if (iq.eq.i_n2) then
	658	DO ilayer=1,nlayer
	659	q(ilayer,iq) = 1
	660	ENDDO
	661	else if (iq.eq.i_ch4_gas) then
	662	ch4ref=0.5 ! default value for vmr ch4
	663	PRINT *,'vmr CH4 (%) ?'
	664	call getin("ch4ref",ch4ref)
	665	write(,) " CH4 ref (%) = ",ch4ref
	666	DO ilayer=1,nlayer
	667	q(ilayer,iq) = 0.01ch4ref16./28.
	668	ENDDO
	669	else if (iq.eq.i_co_gas) then
	670	coref=0.05 ! default value for vmr ch4
	671	PRINT *,'vmr CO (%) ?'
	672	call getin("coref",coref)
	673	write(,) " CO ref (%) = ",coref
	674	DO ilayer=1,nlayer
	675	q(ilayer,iq) = 0.01coref28./28.
	676	ENDDO
[3329]	677	else if ((iq.eq.i_haze10).or.(iq.eq.i_haze30)
	678	& .or.(iq.eq.i_haze).or.(iq.eq.i_haze50)
	679	& .or.(iq.eq.i_haze100)) then
	680	hazeref=0. ! default value for haze mmr
	681	PRINT *,'qhaze (kg/kg) ?'
	682	call getin("hazeref",hazeref)
	683	write(,) " haze ref (kg/kg) = ",hazeref
	684	DO ilayer=1,nlayer
	685	q(ilayer,iq) = hazeref
	686	ENDDO
	687	else if (iq.eq.i_prec_haze) then
	688	prechazeref=0. ! default value for vmr ch4
	689	PRINT *,'qprechaze (kg/kg) ?'
	690	call getin("prechazeref",prechazeref)
	691	write(,) " prec haze ref (kg/kg) = ",prechazeref
	692	DO ilayer=1,nlayer
	693	q(ilayer,iq) = prechazeref
	694	ENDDO
[3258]	695	else
	696	DO ilayer=1,nlayer
	697	q(ilayer,iq) = 0.
	698	ENDDO
	699	endif
[3184]	700	enddo ! of do iq=1,nq
[3258]	701
[3184]	702	endif ! of tracer
	703
	704	c Initialisation pour prendre en compte les vents en 1-D
	705	c ------------------------------------------------------
	706	ptif=2.E+0omegsinlat(1)
	707
[3258]	708
[3184]	709	c vent geostrophique
	710	gru=10. ! default value for gru
	711	PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
	712	call getin("u",gru)
	713	write(,) " u = ",gru
	714	grv=0. !default value for grv
	715	PRINT *,'meridional northward component of the geostrophic',
	716	&' wind (m/s) ?'
	717	call getin("v",grv)
	718	write(,) " v = ",grv
	719
	720	! To be clean, also set vertical winds to zero
	721	w(1:nlayer)=0
	722
	723	c Initialisation des vents au premier pas de temps
	724	DO ilayer=1,nlayer
	725	u(ilayer)=gru
	726	v(ilayer)=grv
	727	ENDDO
	728
	729	c energie cinetique turbulente
	730	DO ilevel=1,nlevel
	731	q2(ilevel)=0.E+0
	732	ENDDO
	733
[3232]	734	c emissivity / surface n2 ice ( + h2o ice??)
[3184]	735	c -------------------------------------------
	736	emis(1)=emissiv ! default value for emissivity
	737	PRINT *,'Emissivity of bare ground ?'
	738	call getin("emis",emis(1))
	739	write(,) " emis = ",emis(1)
[3232]	740	emissiv=emis(1) ! we do this so that condense_n2 sets things to the right
[3184]	741	! value if there is no snow
	742
[3258]	743	if(i_n2.gt.0)then
	744	qsurf(i_n2)=0 ! default value for n2ice
[3232]	745	print*,'Initial n2 ice on the surface (kg.m-2)'
[3258]	746	call getin("n2ice",qsurf(i_n2))
	747	write(,) " n2ice = ",qsurf(i_n2)
	748	IF (qsurf(i_n2).ge.1.E+0) THEN
[3232]	749	! if we have some n2 ice on the surface, change emissivity
[3184]	750	if (latitude(1).ge.0) then ! northern hemisphere
	751	emis(1)=emisice(1)
	752	else ! southern hemisphere
	753	emis(1)=emisice(2)
	754	endif
	755	ENDIF
	756	endif
	757
	758	c calcul des pressions et altitudes en utilisant les niveaux sigma
	759	c ----------------------------------------------------------------
	760
	761	c Vertical Coordinates
	762	c """"""""""""""""""""
	763	hybrid=.true.
	764	PRINT *,'Hybrid coordinates ?'
	765	call getin("hybrid",hybrid)
	766	write(,) " hybrid = ", hybrid
	767
	768
	769	autozlevs=.false.
	770	PRINT *,'Auto-discretise vertical levels ?'
	771	call getin("autozlevs",autozlevs)
	772	write(,) " autozlevs = ", autozlevs
	773
	774	pceil = psurf / 1000.0 ! Pascals
	775	PRINT *,'Ceiling pressure (Pa) ?'
	776	call getin("pceil",pceil)
	777	write(,) " pceil = ", pceil
	778
	779	! Test of incompatibility:
	780	! if autozlevs used, cannot have hybrid too
	781
	782	if (autozlevs.and.hybrid) then
	783	print*,'Cannot use autozlevs and hybrid together!'
	784	call abort
	785	endif
	786
	787	if(autozlevs)then
[3258]	788
[3184]	789	open(91,file="z2sig.def",form='formatted')
	790	read(91,*) Hscale
	791	DO ilayer=1,nlayer-2
	792	read(91,*) Hmax
	793	enddo
	794	close(91)
[3258]	795
	796
[3184]	797	print*,'Hmax = ',Hmax,' km'
	798	print*,'Auto-shifting Hscale to:'
	799	! Hscale = Hmax / log(psurf/100.0)
	800	Hscale = Hmax / log(psurf/pceil)
	801	print*,'Hscale = ',Hscale,' km'
[3258]	802
[3184]	803	! none of this matters if we dont care about zlay
[3258]	804
[3184]	805	endif
	806
	807	call disvert_noterre
	808	! now that disvert has been called, initialize module vertical_layers_mod
	809	call init_vertical_layers(nlayer,preff,scaleheight,
	810	& ap,bp,aps,bps,presnivs,pseudoalt)
	811
	812	if(.not.autozlevs)then
	813	! we want only the scale height from z2sig, in order to compute zlay
	814	open(91,file="z2sig.def",form='formatted')
	815	read(91,*) Hscale
	816	close(91)
	817	endif
	818
	819	! if(autozlevs)then
	820	! open(94,file="Hscale.temp",form='formatted')
	821	! read(94,*) Hscale
	822	! close(94)
	823	! endif
	824
	825	DO ilevel=1,nlevel
	826	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
	827	ENDDO
[3258]	828
[3184]	829	DO ilayer=1,nlayer
	830	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
	831	ENDDO
	832
	833
[3258]	834
[3184]	835	DO ilayer=1,nlayer
	836	! zlay(ilayer)=-300.E+0 rlog(play(ilayer)/plev(1))
	837	! & /g
	838	zlay(ilayer)=-1000.0Hscalelog(play(ilayer)/plev(1))
	839	ENDDO
	840
	841	! endif
	842
	843	c profil de temperature au premier appel
	844	c --------------------------------------
	845	pks=psurf**rcp
	846
	847	c altitude en km dans profile: on divise zlay par 1000
	848	tmp1(0)=0.E+0
	849	DO ilayer=1,nlayer
	850	tmp1(ilayer)=zlay(ilayer)/1000.E+0
	851	ENDDO
	852	call profile(nlayer+1,tmp1,tmp2)
	853
	854	tsurf(1)=tmp2(0)
	855	DO ilayer=1,nlayer
	856	temp(ilayer)=tmp2(ilayer)
	857	ENDDO
	858	print*,"check"
	859	PRINT*,"INPUT SURFACE TEMPERATURE",tsurf(1)
	860	PRINT*,"INPUT TEMPERATURE PROFILE",temp
	861
	862	c Initialisation albedo / inertie du sol
	863	c --------------------------------------
	864	albedodat(1)=0.2 ! default value for albedodat
	865	PRINT *,'Albedo of bare ground ?'
	866	call getin("albedo",albedodat(1))
	867	write(,) " albedo = ",albedodat(1)
	868
	869	inertiedat(1,1)=400 ! default value for inertiedat
	870	PRINT *,'Soil thermal inertia (SI) ?'
	871	call getin("inertia",inertiedat(1,1))
	872	write(,) " inertia = ",inertiedat(1,1)
	873
[3916]	874	c Initialisation n2frac
	875	c --------------------------------------
	876	n2frac(1)=1. ! default value for n2frac
	877
[3184]	878	! Initialize soil properties and temperature
	879	! ------------------------------------------
	880	volcapa=1.e6 ! volumetric heat capacity
[3329]	881	lecttsoil=0 ! default value for lecttsoil
	882	call getin("lecttsoil",lecttsoil)
	883	if (lecttsoil==1) then
	884	OPEN(14,file='proftsoil',status='old',form='formatted',err=101)
	885	DO isoil=1,nsoil
[3672]	886	READ (14,*) tsoil(1,isoil)
[3329]	887	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
	888	ENDDO
	889	GOTO 201
	890	101 STOP'fichier proftsoil inexistant'
	891	201 CONTINUE
	892	CLOSE(14)
	893
[3672]	894	else
[3329]	895	DO isoil=1,nsoil
[3184]	896	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
[3672]	897	tsoil(1,isoil)=tsurf(1) ! soil temperature
[3329]	898	ENDDO
	899	endif
[3184]	900
[3329]	901
[3184]	902	! Initialize depths
	903	! -----------------
	904	do isoil=0,nsoil-1
	905	mlayer(isoil)=3.e-2(2.*(isoil-0.5)) ! mid-layer depth
	906	enddo
	907	do isoil=1,nsoil
	908	layer(isoil)=3.e-2(2.*(isoil-1)) ! layer depth
	909	enddo
	910
[3329]	911	! Initialize haze profile
	912	! ------------------------------------------
	913	if (haze) then
	914
	915	lecthaze=0 ! default value for lecthaze
	916	call getin("lecthaze",lecthaze)
	917	if (lecthaze==1) then
	918	OPEN(15,file='profhaze',status='old',form='formatted',err=301)
	919	DO iq=1,nq
	920	if (iq.eq.i_haze) then
	921	DO ilayer=1,nlayer
	922	READ (15,*) q(ilayer,iq)
	923	ENDDO
	924	endif
	925	ENDDO
	926	GOTO 401
[3406]	927	301 STOP'Problem with profhaze file'
[3329]	928	401 CONTINUE
	929	CLOSE(15)
	930	endif
	931	endif
[3184]	932
[3978]	933	! Initialize gas profiles for clouds
	934	! ------------------------------------------
	935	if (callmufi.and.callmuclouds) then
	936	lectC2H2 = 0 ! default value for lectC2H2
	937	call getin("lectC2H2",lectC2H2)
	938
	939	if (lectC2H2 == 1) then
	940	OPEN(15,file='profC2H2',status='old',form='formatted',err=501)
	941	DO iq = 1, nq
	942	if (iq.eq.i_C2H2_mugas) then
	943	DO ilayer=1,nlayer
	944	READ (15,*) q(ilayer,iq)
	945	ENDDO
	946	endif
	947	ENDDO
	948	GOTO 601
	949	501 STOP'Problem with profC2H2 file'
	950	601 CONTINUE
	951	CLOSE(15)
	952	endif
[3184]	953
[3978]	954	lectC2H6 = 0 ! default value for lectC2H6
	955	call getin("lectC2H6",lectC2H6)
	956
	957	if (lectC2H6 == 1) then
	958	OPEN(15,file='profC2H6',status='old',form='formatted',err=502)
	959	DO iq = 1, nq
	960	if (iq.eq.i_C2H6_mugas) then
	961	DO ilayer=1,nlayer
	962	READ (15,*) q(ilayer,iq)
	963	ENDDO
	964	endif
	965	ENDDO
	966	GOTO 602
	967	502 STOP'Problem with profC2H6 file'
	968	602 CONTINUE
	969	CLOSE(15)
	970	endif
	971
	972	lectC4H2 = 0 ! default value for lectC4H2
	973	call getin("lectC4H2",lectC4H2)
	974
	975	if (lectC4H2 == 1) then
	976	OPEN(15,file='profC4H2',status='old',form='formatted',err=503)
	977	DO iq = 1, nq
	978	if (iq.eq.i_C4H2_mugas) then
	979	DO ilayer=1,nlayer
	980	READ (15,*) q(ilayer,iq)
	981	ENDDO
	982	endif
	983	ENDDO
	984	GOTO 603
	985	503 STOP'Problem with profC4H2 file'
	986	603 CONTINUE
	987	CLOSE(15)
	988	endif
	989
	990	lectC6H6 = 0 ! default value for lectC6H6
	991	call getin("lectC6H6",lectC6H6)
	992
	993	if (lectC6H6 == 1) then
	994	OPEN(15,file='profC6H6',status='old',form='formatted',err=504)
	995	DO iq = 1, nq
	996	if (iq.eq.i_C6H6_mugas) then
	997	DO ilayer=1,nlayer
	998	READ (15,*) q(ilayer,iq)
	999	ENDDO
	1000	endif
	1001	ENDDO
	1002	GOTO 604
	1003	504 STOP'Problem with profC6H6 file'
	1004	604 CONTINUE
	1005	CLOSE(15)
	1006	endif
	1007
	1008	lectHCN = 0 ! default value for lectHCN
	1009	call getin("lectHCN",lectHCN)
	1010
	1011	if (lectHCN == 1) then
	1012	OPEN(15,file='profHCN',status='old',form='formatted',err=505)
	1013	DO iq = 1, nq
	1014	if (iq.eq.i_HCN_mugas) then
	1015	DO ilayer=1,nlayer
	1016	READ (15,*) q(ilayer,iq)
	1017	ENDDO
	1018	endif
	1019	ENDDO
	1020	GOTO 605
	1021	505 STOP'Problem with profHCN file'
	1022	605 CONTINUE
	1023	CLOSE(15)
	1024	endif
	1025	endif ! end of callmufi.and.callmuclouds
	1026
	1027
[3184]	1028	c Write a "startfi" file
	1029	c --------------------
	1030	c This file will be read during the first call to "physiq".
	1031	c It is needed to transfert physics variables to "physiq"...
	1032
	1033	call physdem0("startfi.nc",longitude,latitude,nsoilmx,1,llm,nq,
	1034	& dtphys,real(day0),time,cell_area,
[3541]	1035	& albedodat,zmea,zstd,zsig,zgam,zthe)
[3184]	1036	call physdem1("startfi.nc",nsoilmx,1,llm,nq,
	1037	& dtphys,time,
[3916]	1038	& tsurf,tsoil,inertiedat,emis,albedodat,q2,qsurf,
	1039	& n2frac)
[3184]	1040
	1041	c=======================================================================
[3258]	1042	c BOUCLE TEMPORELLE DU MODELE 1D
[3184]	1043	c=======================================================================
	1044
	1045	firstcall=.true.
	1046	lastcall=.false.
	1047
	1048	DO idt=1,ndt
	1049	IF (idt.eq.ndt) then !test
	1050	lastcall=.true.
	1051	call stellarlong(day*1.0,zls)
	1052	! write(103,) 'Ls=',zls180./pi
	1053	! write(103,) 'Lat=', latitude(1)180./pi
	1054	! write(103,*) 'RunEnd - Atmos. Temp. File'
	1055	! write(103,*) 'RunEnd - Atmos. Temp. File'
	1056	! write(104,) 'Ls=',zls180./pi
	1057	! write(104,*) 'Lat=', latitude(1)
	1058	! write(104,*) 'RunEnd - Atmos. Temp. File'
	1059	ENDIF
	1060
[3258]	1061	c calcul du geopotentiel
[3184]	1062	c ~~~~~~~~~~~~~~~~~~~~~
	1063
	1064
	1065	DO ilayer=1,nlayer
	1066
	1067	! if(autozlevs)then
	1068	! s(ilayer)=(play(ilayer)/psurf)**rcp
	1069	! else
	1070	s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
	1071	! endif
	1072	!s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
	1073	h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
	1074	ENDDO
	1075
	1076	! DO ilayer=1,nlayer
	1077	! s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
	1078	! h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
	1079	! ENDDO
	1080	phi(1)=pksh(1)(1.E+0-s(1))
	1081	DO ilayer=2,nlayer
	1082	phi(ilayer)=phi(ilayer-1)+
	1083	& pks(h(ilayer-1)+h(ilayer)).5E+0
	1084	& *(s(ilayer-1)-s(ilayer))
	1085
	1086	ENDDO
	1087
	1088	c appel de la physique
	1089	c --------------------
	1090
	1091
	1092	CALL physiq (1,llm,nq,
	1093	, firstcall,lastcall,
	1094	, day,time,dtphys,
	1095	, plev,play,phi,
[3258]	1096	, u, v,temp, q,
[3184]	1097	, w,
	1098	C - sorties
	1099	s du, dv, dtemp, dq,dpsurf)
	1100
	1101
	1102	c evolution du vent : modele 1D
	1103	c -----------------------------
[3258]	1104
[3184]	1105	c la physique calcule les derivees temporelles de u et v.
	1106	c on y rajoute betement un effet Coriolis.
	1107	c
	1108	c DO ilayer=1,nlayer
	1109	c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
	1110	c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
	1111	c ENDDO
	1112
	1113	c Pour certain test : pas de coriolis a l'equateur
	1114	c if(latitude(1).eq.0.) then
	1115	DO ilayer=1,nlayer
	1116	du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
	1117	dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
	1118	ENDDO
	1119	c end if
	1120	c
[3258]	1121	c
[3184]	1122	c Calcul du temps au pas de temps suivant
	1123	c ---------------------------------------
	1124	firstcall=.false.
	1125	time=time+dtphys/daysec
	1126	IF (time.gt.1.E+0) then
	1127	time=time-1.E+0
	1128	day=day+1
	1129	ENDIF
	1130
	1131	c calcul des vitesses et temperature au pas de temps suivant
	1132	c ----------------------------------------------------------
	1133
	1134	DO ilayer=1,nlayer
	1135	u(ilayer)=u(ilayer)+dtphys*du(ilayer)
	1136	v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
	1137	if(accelerate_temperature) then
	1138	temp(ilayer)=temp(ilayer)+dtphysdtemp(ilayer)
	1139	& max(1.0,play(ilayer)1e-5)*factor_acceleration
	1140	else
	1141	temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
	1142	endif
	1143	ENDDO
	1144
	1145	c calcul des pressions au pas de temps suivant
	1146	c ----------------------------------------------------------
	1147
	1148	psurf=psurf+dtphys*dpsurf(1) ! evolution de la pression de surface
	1149	DO ilevel=1,nlevel
	1150	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
	1151	ENDDO
	1152	DO ilayer=1,nlayer
	1153	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
	1154	ENDDO
	1155
	1156	c calcul traceur au pas de temps suivant
	1157	c --------------------------------------
	1158
	1159	DO iq = 1, nq
	1160	DO ilayer=1,nlayer
	1161	q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
	1162	ENDDO
	1163	END DO
	1164
	1165	c ========================================================
	1166	c GESTION DES SORTIE
	1167	c ========================================================
	1168	if(saveprofile)then
	1169	OPEN(12,file='profile.out',form='formatted')
	1170	write(12,*) tsurf
	1171	DO ilayer=1,llm
	1172	write(12,*) temp(ilayer) !, play(ilayer) !AS12 only temp so that iprofile=8 can be used
	1173	ENDDO
	1174	CLOSE(12)
	1175	endif
[3329]	1176	! save haze profile
[3672]	1177	if (haze.and.lecthaze.eq.1) then
[3329]	1178	OPEN(16,file='profhaze.out',form='formatted')
	1179	DO iq=1,nq
	1180	if (iq.eq.i_haze) then
	1181	DO ilayer=1,nlayer
	1182	write(16,*) q(ilayer,iq)
	1183	ENDDO
	1184	endif
	1185	ENDDO
	1186	CLOSE(16)
	1187	endif
[3184]	1188
	1189	ENDDO ! fin de la boucle temporelle
	1190
	1191	write(,) "rcm1d: Everything is cool."
	1192
	1193	c ========================================================
	1194	end !rcm1d
	1195
[3258]	1196

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

Download in other formats:

Original Format