Context Navigation

rcm1d.F @ 3811

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

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

Download in other formats:

Original Format