Context Navigation

rcm1d.F @ 3539

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

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

Download in other formats:

Original Format