Context Navigation

rcm1d.F @ 2236

Last change on this file since 2236 was 2116, checked in by jvatant, 6 years ago
Fix an error in tracer management in 1D when running with both mufi and chem. --JVO
Property svn:executable set to ``*
File size: 31.1 KB

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

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

Download in other formats:

Original Format