Context Navigation

rcm1d.F @ 3544

Last change on this file since 3544 was 3542, checked in by tbertrand, 14 months ago
LMDZ.PLUTO Removing nres from planet_mod TB
Property svn:executable set to ``*
File size: 35.7 KB

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

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

Download in other formats:

Original Format