Context Navigation

rcm1d.F @ 724

Last change on this file since 724 was 650, checked in by jleconte, 13 years ago

Corrected the temperature used to differentiate sublimation and evaporation in watersat_grad
Minor name changes in watercommon
Better physical parametrization of the effective radius of liquid and icy water cloud particles in callcorrk

(for radfixed=true)

Added consistency check in inifis
Moved 1d water initialization from physiqu to rcm1d

Property svn:executable set to *

File size: 26.2 KB

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

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

Download in other formats:

Original Format