Context Navigation

rcm1d.F @ 1145

Last change on this file since 1145 was 970, checked in by aslmd, 12 years ago
LMZ.GENERIC rcm1d added a warning about tracer number inconsistency to avoid problems.
Property svn:executable set to ``*
File size: 27.6 KB

Rev	Line
[253]	1	program rcm1d
	2
	3	! to use 'getin'
	4	use ioipsl_getincom
[787]	5	use surfdat_h
	6	use comdiurn_h
	7	use comsaison_h
	8	use comsoil_h
	9	USE comgeomfi_h
[253]	10
	11	implicit none
	12
	13	!==================================================================
	14	!
	15	! Purpose
	16	! -------
	17	! Run the physics package of the universal model in a 1D column.
	18	!
	19	! It can be compiled with a command like (e.g. for 25 layers):
	20	! "makegcm -p std -d 25 rcm1d"
	21	! It requires the files "callphys.def", "z2sig.def",
	22	! "traceur.def", and "run.def" with a line "INCLUDEDEF=callphys.def"
	23	!
	24	! Authors
	25	! -------
	26	! Frederic Hourdin
	27	! R. Fournier
	28	! F. Forget
	29	! F. Montmessin (water ice added)
	30	! R. Wordsworth
[589]	31	! B. Charnay
	32	! A. Spiga
[728]	33	! J. Leconte (2012)
[253]	34	!
	35	!==================================================================
	36
	37	#include "dimensions.h"
[534]	38	#include "paramet.h"
[253]	39	#include "dimphys.h"
	40	#include "callkeys.h"
	41	#include "comcstfi.h"
	42	#include "planete.h"
	43	#include "control.h"
	44	#include "comvert.h"
	45	#include "netcdf.inc"
	46	#include "logic.h"
[534]	47	#include "comgeom.h"
[253]	48
	49	c --------------------------------------------------------------
	50	c Declarations
	51	c --------------------------------------------------------------
	52	c
	53	INTEGER unitstart ! unite d'ecriture de "startfi"
	54	INTEGER nlayer,nlevel,nsoil,ndt
	55	INTEGER ilayer,ilevel,isoil,idt,iq
	56	LOGICAl firstcall,lastcall
	57	c
	58	INTEGER day0 ! date initial (sol ; =0 a Ls=0)
	59	REAL day ! date durant le run
	60	REAL time ! time (0<time<1 ; time=0.5 a midi)
	61	REAL play(nlayermx) ! Pressure at the middle of the layers (Pa)
	62	REAL plev(nlayermx+1) ! intermediate pressure levels (pa)
	63	REAL psurf,tsurf
	64	REAL u(nlayermx),v(nlayermx) ! zonal, meridional wind
	65	REAL gru,grv ! prescribed "geostrophic" background wind
	66	REAL temp(nlayermx) ! temperature at the middle of the layers
	67	REAL q(nlayermx,nqmx) ! tracer mixing ratio (e.g. kg/kg)
	68	REAL qsurf(nqmx) ! tracer surface budget (e.g. kg.m-2)
	69	REAL tsoil(nsoilmx) ! subsurface soik temperature (K)
	70	! REAL co2ice ! co2ice layer (kg.m-2) !not used anymore
	71	integer :: i_co2_ice=0 ! tracer index of co2 ice
[650]	72	integer :: i_h2o_ice=0 ! tracer index of h2o ice
	73	integer :: i_h2o_vap=0 ! tracer index of h2o vapor
[253]	74	REAL emis ! surface layer
	75	REAL q2(nlayermx+1) ! Turbulent Kinetic Energy
	76	REAL zlay(nlayermx) ! altitude estimee dans les couches (km)
	77
	78	c Physical and dynamical tandencies (e.g. m.s-2, K/s, Pa/s)
	79	REAL du(nlayermx),dv(nlayermx),dtemp(nlayermx)
	80	REAL dudyn(nlayermx),dvdyn(nlayermx),dtempdyn(nlayermx)
	81	REAL dpsurf
	82	REAL dq(nlayermx,nqmx)
	83	REAL dqdyn(nlayermx,nqmx)
	84
	85	c Various intermediate variables
	86	! INTEGER thermo
	87	REAL zls
	88	REAL phi(nlayermx),h(nlayermx),s(nlayermx)
	89	REAL pks, ptif, w(nlayermx)
	90	REAL qtotinit, mqtot(nqmx),qtot
	91	INTEGER ierr, aslun
	92	REAL tmp1(0:nlayermx),tmp2(0:nlayermx)
	93	Logical tracerdyn
[787]	94	integer :: nq !=1 ! number of tracers
[253]	95
	96	character*2 str2
	97	character (len=7) :: str7
	98
	99	logical oldcompare, earthhack,saveprofile
	100
	101	! added by RW for zlay computation
	102	real Hscale, Hmax, rho, dz
	103
	104	! added by RW for autozlevs computation
	105	real nu, xx, pMIN, zlev, Htop
	106	real logplevs(nlayermx)
	107
	108	! added by BC
[787]	109	REAL cloudfrac(1,nlayermx)
	110	REAL hice(1),totcloudfrac(1)
[650]	111	REAL qzero1D !initial water amount on the ground
[253]	112
[787]	113	! added by AS to avoid the use of adv trac common
	114	character*20 nametrac(nqmx) ! name of the tracer (no need for adv trac common)
	115
	116
[253]	117	c=======================================================================
	118	c INITIALISATION
	119	c=======================================================================
	120
[787]	121	!! those are defined in surfdat_h.F90
	122	IF (.not. ALLOCATED(albedodat)) ALLOCATE(albedodat(1))
	123	IF (.not. ALLOCATED(phisfi)) ALLOCATE(phisfi(1))
	124	IF (.not. ALLOCATED(zmea)) ALLOCATE(zmea(1))
	125	IF (.not. ALLOCATED(zstd)) ALLOCATE(zstd(1))
	126	IF (.not. ALLOCATED(zsig)) ALLOCATE(zsig(1))
	127	IF (.not. ALLOCATED(zgam)) ALLOCATE(zgam(1))
	128	IF (.not. ALLOCATED(zthe)) ALLOCATE(zthe(1))
	129	IF (.not. ALLOCATED(dryness)) ALLOCATE(dryness(1))
	130	IF (.not. ALLOCATED(watercaptag)) ALLOCATE(watercaptag(1))
	131	!! those are defined in comdiurn_h.F90
	132	IF (.not.ALLOCATED(sinlat)) ALLOCATE(sinlat(1))
	133	IF (.not.ALLOCATED(coslat)) ALLOCATE(coslat(1))
	134	IF (.not.ALLOCATED(sinlon)) ALLOCATE(sinlon(1))
	135	IF (.not.ALLOCATED(coslon)) ALLOCATE(coslon(1))
	136	!! those are defined in comsoil_h.F90
	137	IF (.not.ALLOCATED(layer)) ALLOCATE(layer(nsoilmx))
	138	IF (.not.ALLOCATED(mlayer)) ALLOCATE(mlayer(0:nsoilmx-1))
	139	IF (.not.ALLOCATED(inertiedat)) ALLOCATE(inertiedat(1,nsoilmx))
	140	!! those are defined in comgeomfi_h
	141	IF (.not. ALLOCATED(lati)) ALLOCATE(lati(1))
	142	IF (.not. ALLOCATED(long)) ALLOCATE(long(1))
	143	IF (.not. ALLOCATED(area)) ALLOCATE(area(1))
	144
	145
[526]	146	saveprofile=.false.
[601]	147	saveprofile=.true.
[253]	148
[589]	149	c ----------------------------------------
	150	c Default values (corresponding to Mars)
	151	c ----------------------------------------
[253]	152
	153	pi=2.E+0*asin(1.E+0)
	154
	155	c Parametres Couche limite et Turbulence
	156	c --------------------------------------
	157	z0 = 1.e-2 ! surface roughness (m) ~0.01
	158	emin_turb = 1.e-6 ! energie minimale ~1.e-8
	159	lmixmin = 30 ! longueur de melange ~100
	160
	161	c propriete optiques des calottes et emissivite du sol
	162	c ----------------------------------------------------
	163	emissiv= 0.95 ! Emissivite du sol martien ~.95
	164	emisice(1)=0.95 ! Emissivite calotte nord
	165	emisice(2)=0.95 ! Emissivite calotte sud
	166
	167	albedice(1)=0.5 ! Albedo calotte nord
	168	albedice(2)=0.5 ! Albedo calotte sud
	169	iceradius(1) = 100.e-6 ! mean scat radius of CO2 snow (north)
	170	iceradius(2) = 100.e-6 ! mean scat radius of CO2 snow (south)
	171	dtemisice(1) = 2. ! time scale for snow metamorphism (north)
	172	dtemisice(2) = 2. ! time scale for snow metamorphism (south
	173	hybrid=.false.
	174
	175	c ------------------------------------------------------
[589]	176	c Load parameters from "run.def" and "gases.def"
[253]	177	c ------------------------------------------------------
	178
[589]	179	! check if 'rcm1d.def' file is around
	180	open(90,file='rcm1d.def',status='old',form='formatted',
[253]	181	& iostat=ierr)
	182	if (ierr.ne.0) then
[589]	183	write(,) 'Cannot find required file "rcm1d.def"'
	184	write(,) 'which should contain some input parameters'
[253]	185	write(,) ' ... might as well stop here ...'
	186	stop
	187	else
	188	close(90)
	189	endif
	190
[589]	191	! now, run.def is needed anyway. so we create a dummy temporary one
	192	! for ioipsl to work. if a run.def is already here, stop the
	193	! program and ask the user to do a bit of cleaning
	194	open(90,file='run.def',status='old',form='formatted',
	195	& iostat=ierr)
	196	if (ierr.eq.0) then
	197	close(90)
	198	write(,) 'There is already a run.def file.'
	199	write(,) 'This is not compatible with 1D runs.'
	200	write(,) 'Please remove the file and restart the run.'
	201	write(,) 'Runtime parameters are supposed to be in rcm1d.def'
	202	stop
	203	else
	204	call system('touch run.def')
	205	call system("echo 'INCLUDEDEF=callphys.def' >> run.def")
	206	call system("echo 'INCLUDEDEF=rcm1d.def' >> run.def")
	207	endif
	208
[253]	209	! check if we are going to run with or without tracers
	210	write(,) "Run with or without tracer transport ?"
	211	tracer=.false. ! default value
	212	call getin("tracer",tracer)
	213	write(,) " tracer = ",tracer
	214
[589]	215	! OK. now that run.def has been read once -- any variable is in memory.
	216	! so we can dump the dummy run.def
	217	call system("rm -rf run.def")
	218
[253]	219	! while we're at it, check if there is a 'traceur.def' file
	220	! and preocess it, if necessary. Otherwise initialize tracer names
	221	if (tracer) then
	222	! load tracer names from file 'traceur.def'
	223	open(90,file='traceur.def',status='old',form='formatted',
	224	& iostat=ierr)
	225	if (ierr.eq.0) then
	226	write(,) "rcm1d: Reading file traceur.def"
	227	! read number of tracers:
	228	read(90,*,iostat=ierr) nq
	229	if (ierr.ne.0) then
	230	write(,) "rcm1d: error reading number of tracers"
	231	write(,) " (first line of traceur.def) "
	232	stop
	233	else
	234	! check that the number of tracers is indeed nqmx
	235	if (nq.ne.nqmx) then
	236	write(,) "rcm1d: error, wrong number of tracers:"
	237	write(,) "nq=",nq," whereas nqmx=",nqmx
	238	stop
	239	endif
	240	endif
	241
	242	do iq=1,nq
	243	! minimal version, just read in the tracer names, 1 per line
[787]	244	read(90,*,iostat=ierr) nametrac(iq)
[253]	245	if (ierr.ne.0) then
	246	write(,) 'rcm1d: error reading tracer names...'
	247	stop
	248	endif
	249	enddo !of do iq=1,nq
	250	! check for co2_ice / h2o_ice tracers:
	251	i_co2_ice=0
	252	i_h2o_ice=0
[650]	253	i_h2o_vap=0
[253]	254	do iq=1,nq
[787]	255	if (nametrac(iq)=="co2_ice") then
[253]	256	i_co2_ice=iq
[787]	257	elseif (nametrac(iq)=="h2o_ice") then
[253]	258	i_h2o_ice=iq
[787]	259	elseif (nametrac(iq)=="h2o_vap") then
[650]	260	i_h2o_vap=iq
[253]	261	endif
	262	enddo
	263	!if (i_co2_ice==0) then
	264	! write(,) "rcm1d: error, we need a 'co2_ice' tracer"
	265	! write(,) " (add one to traceur.def)"
	266	! stop
	267	!endif
	268	else
	269	write(,) 'Cannot find required file "traceur.def"'
	270	write(,) ' If you want to run with tracers, I need it'
	271	write(,) ' ... might as well stop here ...'
	272	stop
	273	endif
	274	close(90)
	275
	276	else
[970]	277	nq=nqmx
	278	! Check that tracer boolean is compliant with number of tracers
	279	! -- otherwise there is an error (and more generally we have to be consistent)
	280	if (.not.tracer .and. nq .ge. 1) then
	281	write(,) "------------------------------"
	282	write(,) "rcm1d: You set tracer=.false."
	283	write(,) " But compiled with 1 tracer or more"
	284	write(,) " > If you want tracers, set tracer=.true."
	285	write(,) " > If you do not want tracers, compile with -t 0"
	286	write(,) "------------------------------"
	287	stop
	288	endif
	289	! we still need to set (dummy) tracer names for physdem1
[253]	290	do iq=1,nq
	291	write(str7,'(a1,i2.2)')'q',iq
[787]	292	nametrac(iq)=str7
[253]	293	enddo
	294	! actually, we'll need at least one "co2_ice" tracer
	295	! (for surface CO2 ice)
[787]	296	nametrac(1)="co2_ice"
[253]	297	i_co2_ice=1
	298	endif ! of if (tracer)
	299
[589]	300	!!! We have to check that check_cpp_match is F for 1D computations
	301	!!! We think this check is better than make a particular case for 1D in inifis or calc_cpp_mugaz
	302	check_cpp_match = .false.
	303	call getin("check_cpp_match",check_cpp_match)
	304	if (check_cpp_match) then
	305	print*,"In 1D modeling, check_cpp_match is supposed to be F"
	306	print*,"Please correct callphys.def"
	307	stop
	308	endif
[253]	309
[589]	310	!!! GEOGRAPHICAL INITIALIZATIONS
	311	!!! AREA. useless in 1D
	312	area(1)=1.E+0
	313	aire(1)=area(1) !JL+EM to have access to the area in the diagfi.nc files. area in comgeomfi.h and aire in comgeom.h
[650]	314	!!! GEOPOTENTIAL. useless in 1D because control by surface pressure
[589]	315	phisfi(1)=0.E+0
	316	!!! LATITUDE. can be set.
	317	lati(1)=0 ! default value for lati(1)
	318	PRINT *,'latitude (in degrees) ?'
	319	call getin("latitude",lati(1))
	320	write(,) " latitude = ",lati(1)
	321	lati(1)=lati(1)*pi/180.E+0
	322	!!! LONGITUDE. useless in 1D.
	323	long(1)=0.E+0
	324	long(1)=long(1)*pi/180.E+0
	325
	326	!!! CALL INIFIS
	327	!!! - read callphys.def
	328	!!! - calculate sine and cosine of longitude and latitude
	329	!!! - calculate mugaz and cp
	330	!!! NB: some operations are being done dummily in inifis in 1D
	331	!!! - physical constants: nevermind, things are done allright below
	332	!!! - physical frequency: nevermind, in inifis this is a simple print
[773]	333	cpp=1.d-7 !JL because we divide by cpp in inifis, there may be a more elegant solution
[589]	334	CALL inifis(1,llm,day0,daysec,dtphys,
	335	. lati,long,area,rad,g,r,cpp)
	336	!!! We check everything is OK.
	337	PRINT *,"CHECK"
	338	PRINT *,"--> mugaz = ",mugaz
	339	PRINT *,"--> cpp = ",cpp
	340	r = 8.314511E+0 * 1000.E+0 / mugaz
	341	rcp = r / cpp
	342
	343
	344	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	345	!!!! PLANETARY CONSTANTS !!!!
	346	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	347
	348	g = -99999.
	349	PRINT *,'GRAVITY in m s-2 ?'
	350	call getin("g",g)
	351	IF (g.eq.-99999.) THEN
	352	PRINT *,"STOP. I NEED g IN RUN.DEF."
	353	STOP
	354	ELSE
	355	PRINT *,"--> g = ",g
	356	ENDIF
	357
	358	!rad = -99999.
	359	!PRINT *,'PLANETARY RADIUS in m ?'
	360	!call getin("rad",rad)
	361	!IF (rad.eq.-99999.) THEN
	362	! PRINT *,"STOP. I NEED rad IN RUN.DEF."
	363	! STOP
	364	!ELSE
	365	! PRINT *,"--> rad = ",rad
	366	!ENDIF
	367
	368	daysec = -99999.
	369	PRINT *,'LENGTH OF A DAY in s ?'
	370	call getin("daysec",daysec)
	371	IF (daysec.eq.-99999.) THEN
	372	PRINT *,"STOP. I NEED daysec IN RUN.DEF."
	373	STOP
	374	ELSE
	375	PRINT *,"--> daysec = ",daysec
	376	ENDIF
	377	omeg=4.*asin(1.)/(daysec)
	378	PRINT *,"OK. FROM THIS I WORKED OUT:"
	379	PRINT *,"--> omeg = ",omeg
	380
	381	year_day = -99999.
	382	PRINT *,'LENGTH OF A YEAR in days ?'
	383	call getin("year_day",year_day)
	384	IF (year_day.eq.-99999.) THEN
	385	PRINT *,"STOP. I NEED year_day IN RUN.DEF."
	386	STOP
	387	ELSE
	388	PRINT *,"--> year_day = ",year_day
	389	ENDIF
	390
	391	periastr = -99999.
	392	PRINT *,'MIN DIST STAR-PLANET in AU ?'
	393	call getin("periastr",periastr)
	394	IF (periastr.eq.-99999.) THEN
	395	PRINT *,"STOP. I NEED periastr IN RUN.DEF."
	396	STOP
	397	ELSE
	398	PRINT *,"--> periastr = ",periastr
	399	ENDIF
	400
	401	apoastr = -99999.
[591]	402	PRINT *,'MAX DIST STAR-PLANET in AU ?'
[589]	403	call getin("apoastr",apoastr)
	404	IF (apoastr.eq.-99999.) THEN
	405	PRINT *,"STOP. I NEED apoastr IN RUN.DEF."
	406	STOP
	407	ELSE
	408	PRINT *,"--> apoastr = ",apoastr
	409	ENDIF
	410
	411	peri_day = -99999.
	412	PRINT *,'DATE OF PERIASTRON in days ?'
	413	call getin("peri_day",peri_day)
	414	IF (peri_day.eq.-99999.) THEN
	415	PRINT *,"STOP. I NEED peri_day IN RUN.DEF."
	416	STOP
	417	ELSE IF (peri_day.gt.year_day) THEN
	418	PRINT *,"STOP. peri_day.gt.year_day"
	419	STOP
	420	ELSE
	421	PRINT *,"--> peri_day = ", peri_day
	422	ENDIF
	423
	424	obliquit = -99999.
	425	PRINT *,'OBLIQUITY in deg ?'
	426	call getin("obliquit",obliquit)
	427	IF (obliquit.eq.-99999.) THEN
	428	PRINT *,"STOP. I NEED obliquit IN RUN.DEF."
	429	STOP
	430	ELSE
	431	PRINT *,"--> obliquit = ",obliquit
	432	ENDIF
	433
	434	psurf = -99999.
	435	PRINT *,'SURFACE PRESSURE in Pa ?'
	436	call getin("psurf",psurf)
	437	IF (psurf.eq.-99999.) THEN
	438	PRINT *,"STOP. I NEED psurf IN RUN.DEF."
	439	STOP
	440	ELSE
	441	PRINT *,"--> psurf = ",psurf
	442	ENDIF
	443	!! we need reference pressures
	444	pa=psurf/30.
	445	preff=psurf ! these values are not needed in 1D (are you sure JL12)
	446
	447	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	448	!!!! END PLANETARY CONSTANTS !!!!
	449	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	450
[253]	451	c Date et heure locale du debut du run
	452	c ------------------------------------
	453	c Date (en sols depuis le solstice de printemps) du debut du run
	454	day0 = 0 ! default value for day0
	455	write(,) 'Initial date (in martian sols ; =0 at Ls=0)?'
	456	call getin("day0",day0)
	457	day=float(day0)
	458	write(,) " day0 = ",day0
	459	c Heure de demarrage
	460	time=0 ! default value for time
	461	write(,)'Initial local time (in hours, between 0 and 24)?'
	462	call getin("time",time)
	463	write(,)" time = ",time
	464	time=time/24.E+0 ! convert time (hours) to fraction of sol
	465
[526]	466
[253]	467	c Discretisation (Definition de la grille et des pas de temps)
	468	c --------------
	469	c
	470	nlayer=nlayermx
	471	nlevel=nlayer+1
	472	nsoil=nsoilmx
	473
	474	day_step=48 ! default value for day_step
	475	PRINT *,'Number of time steps per sol ?'
	476	call getin("day_step",day_step)
	477	write(,) " day_step = ",day_step
	478
[534]	479
	480	ecritphy=day_step ! default value for ecritphy
	481	PRINT *,'Nunber of steps between writediagfi ?'
	482	call getin("ecritphy",ecritphy)
	483	write(,) " ecritphy = ",ecritphy
	484
[253]	485	ndt=10 ! default value for ndt
	486	PRINT *,'Number of sols to run ?'
	487	call getin("ndt",ndt)
	488	write(,) " ndt = ",ndt
	489
	490	ndt=ndt*day_step
	491	dtphys=daysec/day_step
[589]	492	write(,)"-------------------------------------"
	493	write(,)"-------------------------------------"
	494	write(,)"--> Physical timestep is ",dtphys
	495	write(,)"-------------------------------------"
	496	write(,)"-------------------------------------"
[253]	497
[526]	498	c output spectrum?
	499	write(,)"Output spectral OLR?"
	500	specOLR=.false.
	501	call getin("specOLR",specOLR)
	502	write(,)" specOLR = ",specOLR
	503
[253]	504	c
	505	c pour le schema d'ondes de gravite
	506	c ---------------------------------
	507	zmea(1)=0.E+0
	508	zstd(1)=0.E+0
	509	zsig(1)=0.E+0
	510	zgam(1)=0.E+0
	511	zthe(1)=0.E+0
	512
	513	c Initialisation des traceurs
	514	c ---------------------------
	515
[787]	516	DO iq=1,nq
[253]	517	DO ilayer=1,nlayer
	518	q(ilayer,iq) = 0.
	519	ENDDO
	520	ENDDO
	521
[787]	522	DO iq=1,nq
[253]	523	qsurf(iq) = 0.
	524	ENDDO
	525
[650]	526	if (water) then
	527	qzero1D = 0.0
	528	qsurf(i_h2o_vap) = qzero1D
	529	endif
	530
[253]	531	c Initialisation pour prendre en compte les vents en 1-D
	532	c ------------------------------------------------------
	533	ptif=2.E+0omegsinlat(1)
	534
	535
	536	c vent geostrophique
	537	gru=10. ! default value for gru
	538	PRINT *,'zonal eastward component of the geostrophic wind (m/s) ?'
	539	call getin("u",gru)
	540	write(,) " u = ",gru
	541	grv=0. !default value for grv
	542	PRINT *,'meridional northward component of the geostrophic',
	543	&' wind (m/s) ?'
	544	call getin("v",grv)
	545	write(,) " v = ",grv
	546
	547	c Initialisation des vents au premier pas de temps
	548	DO ilayer=1,nlayer
	549	u(ilayer)=gru
	550	v(ilayer)=grv
	551	ENDDO
	552
	553	c energie cinetique turbulente
	554	DO ilevel=1,nlevel
	555	q2(ilevel)=0.E+0
	556	ENDDO
	557
	558	c emissivity / surface co2 ice ( + h2o ice??)
	559	c -------------------------------------------
	560	emis=emissiv ! default value for emissivity
	561	PRINT *,'Emissivity of bare ground ?'
	562	call getin("emis",emis)
	563	write(,) " emis = ",emis
	564	emissiv=emis ! we do this so that condense_co2 sets things to the right
	565	! value if there is no snow
	566
	567	if(i_co2_ice.gt.0)then
	568	qsurf(i_co2_ice)=0 ! default value for co2ice
	569	print*,'Initial CO2 ice on the surface (kg.m-2)'
	570	call getin("co2ice",qsurf(i_co2_ice))
	571	write(,) " co2ice = ",qsurf(i_co2_ice)
	572	IF (qsurf(i_co2_ice).ge.1.E+0) THEN
	573	! if we have some CO2 ice on the surface, change emissivity
	574	if (lati(1).ge.0) then ! northern hemisphere
	575	emis=emisice(1)
	576	else ! southern hemisphere
	577	emis=emisice(2)
	578	endif
	579	ENDIF
	580	endif
	581
	582	c calcul des pressions et altitudes en utilisant les niveaux sigma
	583	c ----------------------------------------------------------------
	584
	585	c Vertical Coordinates
	586	c """"""""""""""""""""
	587	hybrid=.true.
	588	PRINT *,'Hybrid coordinates ?'
	589	call getin("hybrid",hybrid)
	590	write(,) " hybrid = ", hybrid
	591
	592
	593	autozlevs=.false.
	594	PRINT *,'Auto-discretise vertical levels ?'
	595	call getin("autozlevs",autozlevs)
	596	write(,) " autozlevs = ", autozlevs
	597
[589]	598	pceil = psurf / 1000.0 ! Pascals
[253]	599	PRINT *,'Ceiling pressure (Pa) ?'
	600	call getin("pceil",pceil)
	601	write(,) " pceil = ", pceil
	602
	603	! Test of incompatibility:
	604	! if autozlevs used, cannot have hybrid too
	605
	606	if (autozlevs.and.hybrid) then
	607	print*,'Cannot use autozlevs and hybrid together!'
	608	call abort
	609	endif
	610
	611	if(autozlevs)then
	612
[305]	613	open(91,file="z2sig.def",form='formatted')
	614	read(91,*) Hscale
	615	DO ilayer=1,nlayer-2
	616	read(91,*) Hmax
	617	enddo
	618	close(91)
[253]	619
	620
	621	print*,'Hmax = ',Hmax,' km'
	622	print*,'Auto-shifting Hscale to:'
	623	! Hscale = Hmax / log(psurf/100.0)
	624	Hscale = Hmax / log(psurf/pceil)
	625	print*,'Hscale = ',Hscale,' km'
	626
	627	! none of this matters if we dont care about zlay
	628
	629	endif
	630
	631	call disvert
	632
	633	if(.not.autozlevs)then
	634	! we want only the scale height from z2sig, in order to compute zlay
	635	open(91,file="z2sig.def",form='formatted')
	636	read(91,*) Hscale
	637	close(91)
	638	endif
	639
	640	! if(autozlevs)then
	641	! open(94,file="Hscale.temp",form='formatted')
	642	! read(94,*) Hscale
	643	! close(94)
	644	! endif
	645
	646	DO ilevel=1,nlevel
	647	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
	648	ENDDO
	649
	650	DO ilayer=1,nlayer
	651	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
	652	ENDDO
	653
	654
	655
	656	DO ilayer=1,nlayer
	657	! zlay(ilayer)=-300.E+0 rlog(play(ilayer)/plev(1))
	658	! & /g
	659	zlay(ilayer)=-1000.0Hscalelog(play(ilayer)/plev(1))
	660	ENDDO
	661
	662	! endif
	663
	664	c profil de temperature au premier appel
	665	c --------------------------------------
	666	pks=psurf**rcp
	667
	668	c altitude en km dans profile: on divise zlay par 1000
	669	tmp1(0)=0.E+0
	670	DO ilayer=1,nlayer
	671	tmp1(ilayer)=zlay(ilayer)/1000.E+0
	672	ENDDO
	673	call profile(nlayer+1,tmp1,tmp2)
	674
	675	tsurf=tmp2(0)
	676	DO ilayer=1,nlayer
	677	temp(ilayer)=tmp2(ilayer)
	678	ENDDO
	679
	680
[589]	681	c Initialisation albedo / inertie du sol
	682	c --------------------------------------
	683	albedodat(1)=0.2 ! default value for albedodat
	684	PRINT *,'Albedo of bare ground ?'
	685	call getin("albedo",albedodat(1))
	686	write(,) " albedo = ",albedodat(1)
[253]	687
[589]	688	inertiedat(1,1)=400 ! default value for inertiedat
	689	PRINT *,'Soil thermal inertia (SI) ?'
	690	call getin("inertia",inertiedat(1,1))
	691	write(,) " inertia = ",inertiedat(1,1)
	692
[253]	693	! Initialize soil properties and temperature
	694	! ------------------------------------------
	695	volcapa=1.e6 ! volumetric heat capacity
	696	DO isoil=1,nsoil
	697	inertiedat(1,isoil)=inertiedat(1,1) ! soil thermal inertia
	698	tsoil(isoil)=tsurf ! soil temperature
	699	ENDDO
	700
	701	! Initialize depths
	702	! -----------------
	703	do isoil=0,nsoil-1
	704	mlayer(isoil)=2.e-4(2.*(isoil-0.5)) ! mid-layer depth
	705	enddo
	706	do isoil=1,nsoil
	707	layer(isoil)=2.e-4(2.*(isoil-1)) ! layer depth
	708	enddo
	709
	710
	711	c Ecriture de "startfi"
	712	c --------------------
	713	c (Ce fichier sera aussitot relu au premier
	714	c appel de "physiq", mais il est necessaire pour passer
	715	c les variables purement physiques a "physiq"...
	716
[787]	717	call physdem1(1,"startfi.nc",long,lati,nsoilmx,nq,
[253]	718	& dtphys,float(day0),
	719	& time,tsurf,tsoil,emis,q2,qsurf,
	720	& area,albedodat,inertiedat,zmea,zstd,zsig,zgam,zthe,
[787]	721	& cloudfrac,totcloudfrac,hice,nametrac)
[253]	722
	723	c=======================================================================
	724	c BOUCLE TEMPORELLE DU MODELE 1D
	725	c=======================================================================
	726
	727	firstcall=.true.
	728	lastcall=.false.
	729
	730	DO idt=1,ndt
[486]	731	IF (idt.eq.ndt) then !test
[253]	732	lastcall=.true.
	733	call stellarlong(day*1.0,zls)
[526]	734	! write(103,) 'Ls=',zls180./pi
	735	! write(103,) 'Lat=', lati(1)180./pi
	736	! write(103,*) 'RunEnd - Atmos. Temp. File'
	737	! write(103,*) 'RunEnd - Atmos. Temp. File'
	738	! write(104,) 'Ls=',zls180./pi
	739	! write(104,*) 'Lat=', lati(1)
	740	! write(104,*) 'RunEnd - Atmos. Temp. File'
[253]	741	ENDIF
	742
	743	c calcul du geopotentiel
	744	c ~~~~~~~~~~~~~~~~~~~~~
	745
	746
[586]	747	DO ilayer=1,nlayer
[253]	748
	749	! if(autozlevs)then
	750	! s(ilayer)=(play(ilayer)/psurf)**rcp
	751	! else
[586]	752	s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
[253]	753	! endif
	754	!s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
[586]	755	h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
	756	ENDDO
	757
[253]	758	! DO ilayer=1,nlayer
	759	! s(ilayer)=(aps(ilayer)/psurf+bps(ilayer))**rcp
	760	! h(ilayer)=cpptemp(ilayer)/(pkss(ilayer))
	761	! ENDDO
	762	phi(1)=pksh(1)(1.E+0-s(1))
	763	DO ilayer=2,nlayer
	764	phi(ilayer)=phi(ilayer-1)+
	765	& pks(h(ilayer-1)+h(ilayer)).5E+0
	766	& *(s(ilayer-1)-s(ilayer))
	767
	768	ENDDO
	769
	770	c appel de la physique
	771	c --------------------
	772
	773
[787]	774	CALL physiq (1,llm,nq,
	775	. nametrac,
[253]	776	, firstcall,lastcall,
	777	, day,time,dtphys,
	778	, plev,play,phi,
	779	, u, v,temp, q,
	780	, w,
	781	C - sorties
	782	s du, dv, dtemp, dq,dpsurf,tracerdyn)
	783
	784
	785	c evolution du vent : modele 1D
	786	c -----------------------------
	787
	788	c la physique calcule les derivees temporelles de u et v.
	789	c on y rajoute betement un effet Coriolis.
	790	c
	791	c DO ilayer=1,nlayer
	792	c du(ilayer)=du(ilayer)+ptif*(v(ilayer)-grv)
	793	c dv(ilayer)=dv(ilayer)+ptif*(-u(ilayer)+gru)
	794	c ENDDO
	795
	796	c Pour certain test : pas de coriolis a l'equateur
	797	c if(lati(1).eq.0.) then
	798	DO ilayer=1,nlayer
	799	du(ilayer)=du(ilayer)+ (gru-u(ilayer))/1.e4
	800	dv(ilayer)=dv(ilayer)+ (grv-v(ilayer))/1.e4
	801	ENDDO
	802	c end if
	803	c
	804	c
	805	c Calcul du temps au pas de temps suivant
	806	c ---------------------------------------
	807	firstcall=.false.
	808	time=time+dtphys/daysec
	809	IF (time.gt.1.E+0) then
	810	time=time-1.E+0
	811	day=day+1
	812	ENDIF
	813
	814	c calcul des vitesses et temperature au pas de temps suivant
	815	c ----------------------------------------------------------
	816
	817	DO ilayer=1,nlayer
	818	u(ilayer)=u(ilayer)+dtphys*du(ilayer)
	819	v(ilayer)=v(ilayer)+dtphys*dv(ilayer)
	820	temp(ilayer)=temp(ilayer)+dtphys*dtemp(ilayer)
	821	ENDDO
	822
	823	c calcul des pressions au pas de temps suivant
	824	c ----------------------------------------------------------
	825
	826	psurf=psurf+dtphys*dpsurf ! evolution de la pression de surface
	827	DO ilevel=1,nlevel
	828	plev(ilevel)=ap(ilevel)+psurf*bp(ilevel)
	829	ENDDO
	830	DO ilayer=1,nlayer
	831	play(ilayer)=aps(ilayer)+psurf*bps(ilayer)
	832	ENDDO
	833
	834	c calcul traceur au pas de temps suivant
	835	c --------------------------------------
	836
[787]	837	DO iq = 1, nq
[253]	838	DO ilayer=1,nlayer
	839	q(ilayer,iq)=q(ilayer,iq)+dtphys*dq(ilayer,iq)
	840	ENDDO
	841	END DO
	842
	843	c ========================================================
	844	c GESTION DES SORTIE
	845	c ========================================================
	846	if(saveprofile)then
	847	OPEN(12,file='profile.out',form='formatted')
[601]	848	write(12,*) tsurf
[253]	849	DO ilayer=1,nlayermx
[601]	850	write(12,*) temp(ilayer) !, play(ilayer) !AS12 only temp so that iprofile=8 can be used
[253]	851	ENDDO
	852	CLOSE(12)
	853	endif
	854
	855
[601]	856	ENDDO ! fin de la boucle temporelle
	857
[253]	858	c ========================================================
	859	end !rcm1d
	860
	861	c***********************************************************************
	862	c***********************************************************************
	863	c Subroutines Bidons utilise seulement en 3D, mais
	864	c necessaire a la compilation de rcm1d en 1D
	865
	866	subroutine gr_fi_dyn
	867	RETURN
	868	END
	869
	870	c***********************************************************************
	871	c***********************************************************************
	872
	873	#include "../dyn3d/disvert.F"

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

Download in other formats:

Original Format