Context Navigation

rcm1d.F @ 1533

Last change on this file since 1533 was 1533, checked in by mturbet, 9 years ago
Tracers profiles now available for input
Property svn:executable set to ``*
File size: 31.6 KB

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

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

Download in other formats:

Original Format