Context Navigation

rcm1d.F @ 1802

Last change on this file since 1802 was 1802, checked in by bclmd, 7 years ago
Adding photochemistry to LMDZ Generic
Property svn:executable set to ``*
File size: 33.1 KB

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

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

Download in other formats:

Original Format