Context Navigation

rcm1d.F @ 2324

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

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

Download in other formats:

Original Format