Context Navigation

rcm1d.F @ 2613

Last change on this file since 2613 was 2542, checked in by aslmd, 3 years ago

Generic GCM:

Large update of the chemical modules

Read chemical network from input files
Init chemistry with ModernTrac?
Photolysis online calculation

Property svn:executable set to *

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

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

Download in other formats:

Original Format