Context Navigation

inichim_newstart.F90 @ 1862

Last change on this file since 1862 was 1806, checked in by bclmd, 7 years ago
inichim_newstart with Mars case by default
Property svn:executable set to ``*
File size: 20.8 KB

Rev	Line
[1796]	1	subroutine inichim_newstart(ngrid, nq, pq, qsurf, ps, &
	2	flagh2o, flagthermo)
	3
	4	use tracer_h
	5	USE comvert_mod, ONLY: aps,bps
	6	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev
	7	use callkeys_mod
	8	use datafile_mod
	9
	10	implicit none
	11
	12	!=======================================================================
	13	!
	14	! Purpose:
	15	! --------
	16	!
	17	! Initialization of the chemistry for use in the building of a new start file
	18	! used by program newstart.F
	19	! also used by program testphys1d.F
	20	!
	21	! Authors:
	22	! -------
	23	! Initial version 11/2002 by Sebastien Lebonnois
	24	! Revised 07/2003 by Monica Angelats-i-Coll to use input files
	25	! Modified 10/2008 Identify tracers by their names Ehouarn Millour
	26	! Modified 11/2011 Addition of methane Franck Lefevre
	27	! Rewritten 04/2012 Franck Lefevre
	28	!
	29	! Arguments:
	30	! ----------
	31	!
	32	! pq(nbp_lon+1,nbp_lat,nbp_lev,nq) Advected fields, ie chemical species here
	33	! qsurf(ngrid,nq) Amount of tracer on the surface (kg/m2)
	34	! ps(nbp_lon+1,nbp_lat) Surface pressure (Pa)
	35	! flagh2o flag for initialisation of h2o (1: yes / 0: no)
	36	! flagthermo flag for initialisation of thermosphere only (1: yes / 0: no)
	37	!
	38	!=======================================================================
	39
	40
	41	! inputs :
	42
	43	integer,intent(in) :: ngrid ! number of atmospheric columns in the physics
	44	integer,intent(in) :: nq ! number of tracers
	45	real,intent(in) :: ps(nbp_lon+1,nbp_lat) ! surface pressure in the gcm (Pa)
	46	integer,intent(in) :: flagh2o ! flag for h2o initialisation
	47	integer,intent(in) :: flagthermo ! flag for thermosphere initialisation only
	48
	49	! outputs :
	50
	51	real,intent(out) :: pq(nbp_lon+1,nbp_lat,nbp_lev,nq) ! advected fields, ie chemical species
	52	real,intent(out) :: qsurf(ngrid,nq) ! surface values (kg/m2) of tracers
	53
	54	! local :
	55
	56	integer :: iq, i, j, l, n, nbqchem
	57	integer :: count, ierr, dummy
	58	real :: mmean(nbp_lon+1,nbp_lat,nbp_lev) ! mean molecular mass (g)
	59	real :: pgcm ! pressure at each layer in the gcm (Pa)
	60
	61	integer, parameter :: nalt = 252 ! number of altitudes in the initialization files
	62	integer :: nspe ! number of species in the initialization files
	63	integer, allocatable :: niq(:) ! local index of species in initialization files
	64	real, dimension(nalt) :: tinit, zzfile ! temperature in initialization files
	65	real, dimension(nalt) :: pinit ! pressure in initialization files
	66	real, dimension(nalt) :: densinit ! total number density in initialization files
	67	real, allocatable :: vmrinit(:,:) ! mixing ratios in initialization files
	68	real, allocatable :: vmrint(:) ! mixing ratio interpolated onto the gcm vertical grid
	69	real :: vmr
	70
	71	character(len=20) :: txt ! to store some text
	72	logical :: flagnitro ! checks if N species present
	73
	74	! 1. identify tracers by their names: (and set corresponding values of mmol)
	75
	76	! 1.1 initialize tracer indexes to zero:
	77	! nqmx=nq ! initialize value of nqmx
	78
	79	do iq = 1,nq
	80	igcm_dustbin(iq) = 0
	81	end do
	82
	83	! igcm_dust_mass = 0
	84	! igcm_dust_number = 0
	85	! igcm_ccn_mass = 0
	86	! igcm_ccn_number = 0
	87	igcm_h2o_vap = 0
	88	igcm_h2o_ice = 0
	89	igcm_co2 = 0
	90	igcm_co = 0
	91	igcm_o = 0
	92	igcm_o1d = 0
	93	igcm_o2 = 0
	94	igcm_o3 = 0
	95	igcm_h = 0
	96	igcm_h2 = 0
	97	igcm_oh = 0
	98	igcm_ho2 = 0
	99	igcm_h2o2 = 0
	100	igcm_ch4 = 0
	101	igcm_n2 = 0
	102	igcm_ar = 0
	103	igcm_ar_n2 = 0
	104	igcm_ch3 = 0
	105	igcm_ch = 0
	106	igcm_3ch2 = 0
	107	igcm_1ch2 = 0
	108	igcm_cho = 0
	109	igcm_ch2o = 0
	110	igcm_c = 0
	111	igcm_c2 = 0
	112	igcm_c2h = 0
	113	igcm_c2h2 = 0
	114	igcm_c2h3 = 0
	115	igcm_c2h4 = 0
	116	igcm_c2h6 = 0
	117	igcm_ch2co = 0
	118	igcm_ch3co = 0
	119	igcm_hcaer = 0
	120
	121	! 1.2 find dust tracers
	122	count = 0
	123	!
	124	! if (dustbin > 0) then
	125	! do iq = 1,nq
	126	! txt = " "
	127	! write(txt,'(a4,i2.2)') 'dust', count + 1
	128	! if (noms(iq) == txt) then
	129	! count = count + 1
	130	! igcm_dustbin(count) = iq
	131	! mmol(iq) = 100.
	132	! end if
	133	! end do !do iq=1,nq
	134	! end if ! of if (dustbin.gt.0)
	135	!
	136	! if (doubleq) then
	137	! do iq = 1,nq
	138	! if (noms(iq) == "dust_mass") then
	139	! igcm_dust_mass = iq
	140	! count = count + 1
	141	! end if
	142	! if (noms(iq) == "dust_number") then
	143	! igcm_dust_number = iq
	144	! count = count + 1
	145	! end if
	146	! end do
	147	! end if ! of if (doubleq)
	148	!
	149	! if (scavenging) then
	150	! do iq = 1,nq
	151	! if (noms(iq) == "ccn_mass") then
	152	! igcm_ccn_mass = iq
	153	! count = count + 1
	154	! end if
	155	! if (noms(iq) == "ccn_number") then
	156	! igcm_ccn_number = iq
	157	! count = count + 1
	158	! end if
	159	! end do
	160	! end if ! of if (scavenging)
	161	!
	162	! if (submicron) then
	163	! do iq=1,nq
	164	! if (noms(iq) == "dust_submicron") then
	165	! igcm_dust_submicron = iq
	166	! mmol(iq) = 100.
	167	! count = count + 1
	168	! end if
	169	! end do
	170	! end if ! of if (submicron)
	171
	172	! 1.3 find chemistry and water tracers
	173
	174	nbqchem = 0
	175
	176	do iq = 1,nq
	177	if (noms(iq) == "co2") then
	178	igcm_co2 = iq
	179	mmol(igcm_co2) = 44.
	180	count = count + 1
	181	nbqchem = nbqchem + 1
	182	end if
	183	if (noms(iq) == "co") then
	184	igcm_co = iq
	185	mmol(igcm_co) = 28.
	186	count = count + 1
	187	nbqchem = nbqchem + 1
	188	end if
	189	if (noms(iq) == "o") then
	190	igcm_o = iq
	191	mmol(igcm_o) = 16.
	192	count = count + 1
	193	nbqchem = nbqchem + 1
	194	end if
	195	if (noms(iq) == "o1d") then
	196	igcm_o1d = iq
	197	mmol(igcm_o1d) = 16.
	198	count = count + 1
	199	nbqchem = nbqchem + 1
	200	end if
	201	if (noms(iq) == "o2") then
	202	igcm_o2 = iq
	203	mmol(igcm_o2) = 32.
	204	count = count + 1
	205	nbqchem = nbqchem + 1
	206	end if
	207	if (noms(iq) == "o3") then
	208	igcm_o3 = iq
	209	mmol(igcm_o3) = 48.
	210	count = count + 1
	211	nbqchem = nbqchem + 1
	212	end if
	213	if (noms(iq) == "h") then
	214	igcm_h = iq
	215	mmol(igcm_h) = 1.
	216	count = count + 1
	217	nbqchem = nbqchem + 1
	218	end if
	219	if (noms(iq) == "h2") then
	220	igcm_h2 = iq
	221	mmol(igcm_h2) = 2.
	222	count = count + 1
	223	nbqchem = nbqchem + 1
	224	end if
	225	if (noms(iq) == "oh") then
	226	igcm_oh = iq
	227	mmol(igcm_oh) = 17.
	228	count = count + 1
	229	nbqchem = nbqchem + 1
	230	end if
	231	if (noms(iq) == "ho2") then
	232	igcm_ho2 = iq
	233	mmol(igcm_ho2) = 33.
	234	count = count + 1
	235	nbqchem = nbqchem + 1
	236	end if
	237	if (noms(iq) == "h2o2") then
	238	igcm_h2o2 = iq
	239	mmol(igcm_h2o2) = 34.
	240	count = count + 1
	241	nbqchem = nbqchem + 1
	242	end if
	243	if (noms(iq) == "ch4") then
	244	igcm_ch4 = iq
	245	mmol(igcm_ch4) = 16.
	246	count = count + 1
	247	nbqchem = nbqchem + 1
	248	end if
	249	if (noms(iq) == "n2") then
	250	igcm_n2 = iq
	251	mmol(igcm_n2) = 28.
	252	count = count + 1
	253	nbqchem = nbqchem + 1
	254	end if
	255	if (noms(iq) == "n") then
	256	igcm_n = iq
	257	mmol(igcm_n) = 14.
	258	count = count + 1
	259	nbqchem = nbqchem + 1
	260	end if
	261	if (noms(iq) == "n2d") then
	262	igcm_n2d = iq
	263	mmol(igcm_n2d) = 14.
	264	count = count + 1
	265	nbqchem = nbqchem + 1
	266	end if
	267	if (noms(iq) == "no") then
	268	igcm_no = iq
	269	mmol(igcm_no) = 30.
	270	count = count + 1
	271	nbqchem = nbqchem + 1
	272	end if
	273	if (noms(iq) == "no2") then
	274	igcm_no2 = iq
	275	mmol(igcm_no2) = 46.
	276	count = count + 1
	277	nbqchem = nbqchem + 1
	278	end if
	279	if (noms(iq) == "h2o_vap") then
	280	igcm_h2o_vap = iq
	281	mmol(igcm_h2o_vap) = 18.
	282	count = count + 1
	283	nbqchem = nbqchem + 1
	284	end if
	285	if (noms(iq) == "h2o_ice") then
	286	igcm_h2o_ice = iq
	287	mmol(igcm_h2o_ice) = 18.
	288	count = count + 1
	289	nbqchem = nbqchem + 1
	290	end if
	291
	292	if (noms(iq).eq."ch3") then
	293	igcm_ch3=iq
	294	mmol(igcm_ch3)=15.
	295	count=count+1
	296	nbqchem = nbqchem + 1
	297	endif
	298	if (noms(iq).eq."ch") then
	299	igcm_ch=iq
	300	mmol(igcm_ch)=13.
	301	count=count+1
	302	nbqchem = nbqchem + 1
	303	endif
	304	if (noms(iq).eq."3ch2") then
	305	igcm_3ch2=iq
	306	mmol(igcm_3ch2)=14.
	307	count=count+1
	308	nbqchem = nbqchem + 1
	309	endif
	310	if (noms(iq).eq."1ch2") then
	311	igcm_1ch2=iq
	312	mmol(igcm_1ch2)=14.
	313	count=count+1
	314	nbqchem = nbqchem + 1
	315	endif
	316	if (noms(iq).eq."cho") then
	317	igcm_cho=iq
	318	mmol(igcm_cho)=29.
	319	count=count+1
	320	nbqchem = nbqchem + 1
	321	endif
	322	if (noms(iq).eq."ch2o") then
	323	igcm_ch2o=iq
	324	mmol(igcm_ch2o)=30.
	325	count=count+1
	326	nbqchem = nbqchem + 1
	327	endif
	328	if (noms(iq).eq."ch3o") then
	329	igcm_ch3o=iq
	330	mmol(igcm_ch3o)=31.
	331	count=count+1
	332	nbqchem = nbqchem + 1
	333	endif
	334	if (noms(iq).eq."c") then
	335	igcm_c=iq
	336	mmol(igcm_c)=12.
	337	count=count+1
	338	nbqchem = nbqchem + 1
	339	endif
	340	if (noms(iq).eq."c2") then
	341	igcm_c2=iq
	342	mmol(igcm_c2)=24.
	343	count=count+1
	344	nbqchem = nbqchem + 1
	345	endif
	346	if (noms(iq).eq."c2h") then
	347	igcm_c2h=iq
	348	mmol(igcm_c2h)=25.
	349	count=count+1
	350	nbqchem = nbqchem + 1
	351	endif
	352	if (noms(iq).eq."c2h2") then
	353	igcm_c2h2=iq
	354	mmol(igcm_c2h2)=26.
	355	count=count+1
	356	nbqchem = nbqchem + 1
	357	endif
	358	if (noms(iq).eq."c2h3") then
	359	igcm_c2h3=iq
	360	mmol(igcm_c2h3)=27.
	361	count=count+1
	362	nbqchem = nbqchem + 1
	363	endif
	364	if (noms(iq).eq."c2h4") then
	365	igcm_c2h4=iq
	366	mmol(igcm_c2h4)=28.
	367	count=count+1
	368	nbqchem = nbqchem + 1
	369	endif
	370	if (noms(iq).eq."c2h6") then
	371	igcm_c2h6=iq
	372	mmol(igcm_c2h6)=30.
	373	count=count+1
	374	nbqchem = nbqchem + 1
	375	endif
	376	if (noms(iq).eq."ch2co") then
	377	igcm_ch2co=iq
	378	mmol(igcm_ch2co)=42.
	379	count=count+1
	380	nbqchem = nbqchem + 1
	381	endif
	382	if (noms(iq).eq."ch3co") then
	383	igcm_ch3co=iq
	384	mmol(igcm_ch3co)=43.
	385	count=count+1
	386	nbqchem = nbqchem + 1
	387	endif
	388	if (noms(iq).eq."hcaer") then
	389	igcm_hcaer=iq
	390	mmol(igcm_hcaer)=50.
	391	count=count+1
	392	nbqchem = nbqchem + 1
	393	endif
	394	if (noms(iq) == "ar") then
	395	igcm_ar = iq
	396	mmol(igcm_ar) = 40.
	397	count = count + 1
	398	nbqchem = nbqchem + 1
	399	end if
	400
	401
	402
	403	! 1.5 find idealized non-condensible tracer
	404
	405	if (noms(iq) == "Ar_N2") then
	406	igcm_ar_n2 = iq
	407	mmol(igcm_ar_n2) = 30.
	408	count = count + 1
	409	end if
	410
	411	end do ! of do iq=1,nq
	412
	413	! 1.6 check that we identified all tracers:
	414
	415	if (count /= nq) then
	416	write(,) "inichim_newstart: found only ",count," tracers"
	417	write(,) " expected ",nq
	418	do iq = 1,count
	419	write(,) ' ', iq, ' ', trim(noms(iq))
	420	end do
	421	stop
	422	else
	423	write(,) "inichim_newstart: found all expected tracers"
	424	do iq = 1,nq
	425	write(,) ' ', iq, ' ', trim(noms(iq))
	426	end do
	427	end if
	428
	429	! 1.7 check if nitrogen species are present:
	430
	431	if(igcm_no == 0) then
	432	!check that no N species is in traceur.def
	433	if(igcm_n /= 0 .or. igcm_no2 /= 0 .or. igcm_n2d /= 0) then
	434	write(,)'inichim_newstart error:'
	435	write(,)'N, NO2 and/or N2D are in traceur.def, but not NO'
	436	write(,)'stop'
	437	stop
	438	endif
	439	flagnitro = .false.
	440	nspe = 14
	441	else
	442	!check that all N species are in traceur.def
	443	if(igcm_n == 0 .or. igcm_no2 == 0 .or. igcm_n2d == 0) then
	444	write(,)'inichim_newstart error:'
	445	write(,)'if NO is in traceur.def, N, NO2 and N2D must also be'
	446	write(,)'stop'
	447	stop
	448	endif
	449	flagnitro = .true.
[1806]	450	nspe = 18
[1796]	451	endif
	452
	453	! 1.8 allocate arrays
	454
	455	allocate(niq(nspe))
	456	allocate(vmrinit(nalt,nspe))
	457	allocate(vmrint(nspe))
	458
	459	! 2. load in chemistry data for initialization:
	460
	461	! order of major species in initialization file:
	462	!
	463	! 1: co2
	464	! 2: ar
	465	! 3: n2
	466	! 4: o2
	467	! 5: co
	468	! 6: o
	469	! 7: h2
	470	!
	471	! order of minor species in initialization file:
	472	!
	473	! 1: h
	474	! 2: oh
	475	! 3: ho2
	476	! 4: h2o
	477	! 5: h2o2
	478	! 6: o1d
	479	! 7: o3
	480	!
	481	! order of nitrogen species in initialization file:
	482	!
	483	! 1: n
	484	! 2: no
	485	! 3: no2
	486	! 4: n2d
	487
	488	! major species:
	489
	490	niq(1) = igcm_co2
[1806]	491	niq(2) = igcm_ar
	492	niq(3) = igcm_n2
	493	niq(4) = igcm_o2
	494	niq(5) = igcm_co
	495	niq(6) = igcm_o
	496	niq(7) = igcm_h2
[1796]	497
	498	! minor species:
	499
[1806]	500	niq(8) = igcm_h
	501	niq(9) = igcm_oh
	502	niq(10) = igcm_ho2
	503	niq(11) = igcm_h2o_vap
	504	niq(12) = igcm_h2o2
	505	niq(13) = igcm_o1d
	506	niq(14) = igcm_o3
[1796]	507
	508	! nitrogen species:
	509
	510	if (flagnitro) then
[1806]	511	niq(15) = igcm_n
	512	niq(16) = igcm_no
	513	niq(17) = igcm_no2
	514	niq(18) = igcm_n2d
[1796]	515	end if
	516
	517
	518
	519
	520	! carbon species:
[1806]	521	! niq(18) = igcm_ch4
	522	! niq(19) = igcm_ch3
	523	! niq(20) = igcm_ch
	524	! niq(21) = igcm_1ch2
	525	! niq(22) = igcm_3ch2
	526	! niq(23) = igcm_cho
	527	! niq(24) = igcm_ch2o
	528	! niq(25) = igcm_ch3o
	529	! niq(26) = igcm_c
	530	! niq(27) = igcm_c2
	531	! niq(28) = igcm_c2h
	532	! niq(29) = igcm_c2h2
	533	! niq(30) = igcm_c2h3
	534	! niq(31) = igcm_c2h4
	535	! niq(32) = igcm_c2h6
	536	! niq(33) = igcm_ch2co
	537	! niq(34) = igcm_ch3co
	538	! niq(35) = igcm_hcaer
[1796]	539
	540
	541	! 2.1 open initialization files
	542	open(210, iostat=ierr,file=trim(datadir)//'/atmosfera_LMD_may.dat')
	543	if (ierr /= 0) then
	544	write(,)'Error : cannot open file atmosfera_LMD_may.dat '
	545	write(,)'(in aeronomars/inichim_newstart.F)'
	546	write(,)'It should be in :', trim(datadir),'/'
	547	write(,)'1) You can change this path in callphys.def with'
	548	write(,)' datadir=/path/to/datafiles/'
	549	write(,)'2) If necessary atmosfera_LMD_may.dat (and others)'
	550	write(,)' can be obtained online on:'
	551	write(,)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
	552	stop
	553	end if
	554	open(220, iostat=ierr,file=trim(datadir)//'/atmosfera_LMD_min.dat')
	555	if (ierr /= 0) then
	556	write(,)'Error : cannot open file atmosfera_LMD_min.dat '
	557	write(,)'(in aeronomars/inichim_newstart.F)'
	558	write(,)'It should be in :', trim(datadir),'/'
	559	write(,)'1) You can change this path in callphys.def with'
	560	write(,)' datadir=/path/to/datafiles/'
	561	write(,)'2) If necessary atmosfera_LMD_min.dat (and others)'
	562	write(,)' can be obtained online on:'
	563	write(,)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
	564	stop
	565	end if
	566	if(flagnitro) then
	567	open(230, iostat=ierr,file=trim(datadir)//'/atmosfera_LMD_nitr.dat')
	568	if (ierr.ne.0) then
	569	write(,)'Error : cannot open file atmosfera_LMD_nitr.dat '
	570	write(,)'(in aeronomars/inichim_newstart.F)'
	571	write(,)'It should be in :', datadir
	572	write(,)'1) You can change this directory address in '
	573	write(,)' file phymars/datafile.h'
	574	write(,)'2) If necessary atmosfera_LMD_nitr.dat (and others)'
	575	write(,)' can be obtained online on:'
	576	write(,)' http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
	577	STOP
	578	endif
	579	endif ! Of if(flagnitro)
	580
	581	! 2.2 read initialization files
	582
	583	! major species
	584
	585	read(210,*)
	586	do l = 1,nalt
	587	read(210,*) dummy, tinit(l), pinit(l), densinit(l), &
	588	(vmrinit(l,n), n = 1,7)
	589	pinit(l) = pinit(l)*100. ! conversion in Pa
	590	pinit(l) = log(pinit(l)) ! for the vertical interpolation
	591	end do
	592	close(210)
	593
	594	! minor species
	595
	596	read(220,*)
	597	do l = 1,nalt
	598	read(220,*) dummy, (vmrinit(l,n), n = 8,14)
	599	end do
	600	close(220)
	601
	602	! nitrogen species
	603
	604	if (flagnitro) then
	605	read(230,*)
	606	do l = 1,nalt
	607	read(230,*) dummy, (vmrinit(l,n), n = 15,18)
	608	end do
	609	close(230)
	610	end if
	611
[1806]	612	! initialization for the early eath
	613	if (1.eq.0) then
[1796]	614	do l = 1,nalt
	615	vmrinit(l,:)=0.0
	616	vmrinit(l,1)=0.01 !co2
	617	vmrinit(l,2)=0.989 !n2
	618	! vmrinit(l,3)=2.e-17/mmol(niq(3))*28 !o2
	619	! vmrinit(l,4)=3.8e-6/mmol(niq(4))*28 !co
	620	! vmrinit(l,5)=4.e-14/mmol(niq(5))*28 !o
	621	! vmrinit(l,6)=1.3e-7/mmol(niq(6))*28 !h2
	622	vmrinit(l,6)=1e-3
	623	! vmrinit(l,7)=5.e-16/mmol(niq(7))*28 !h
	624	! vmrinit(l,8)=2.e-17/mmol(niq(8))*28 !oh
	625	! vmrinit(l,9)=1.e-17/mmol(niq(9))*28 !ho2
	626	vmrinit(l,10)=1e-6 !h2o
	627	! vmrinit(l,11)=2.e-20/mmol(niq(11))*28 !h2o2
	628	! vmrinit(l,12)=0. !o1d
	629	! vmrinit(l,13)=3.e-22/mmol(niq(13))*28 !o3
	630
	631
	632	vmrinit(l,18)=1.0e-3 !ch4
	633	! vmrinit(l,19)=1.3e-12/mmol(niq(19))*28 !ch3
	634	! vmrinit(l,23)=1.e-12/mmol(niq(23))*28 !cho
	635	! vmrinit(l,24)=2.7e-11/mmol(niq(24))*28 !ch2o
	636	! vmrinit(l,25)=2.e-9/mmol(niq(25))*28 !ch3o
	637	! vmrinit(l,32)=2.e-7/mmol(niq(32))*28 !c2h6
	638	! vmrinit(l,33)=5.e-12/mmol(niq(33))*28 !ch2co
	639	! vmrinit(l,34)=1.e-13/mmol(niq(34))*28 !ch3co
	640
	641
	642
	643	! pinit(l)=aps(l) + bps(l)*ps
	644	! vmrinit(l,18)=2e-3*min(pinit(l)/100.,1.) ! decrease with scale
	645	! height above 1 hpa
	646	vmrinit(l,2)=0.0
	647	vmrinit(l,2)=1-sum(vmrinit(l,:)) !n2
	648	! vmrinit(l,4)=0.1
	649	! vmrinit(l,7)=0.001
	650	end do
	651	endif
	652
	653
	654	! 3. initialization of tracers
	655
	656	do i = 1,nbp_lon+1
	657	do j = 1,nbp_lat
	658	do l = 1,nbp_lev
	659
	660	pgcm = aps(l) + bps(l)*ps(i,j) ! gcm pressure
	661	pgcm = log(pgcm) ! for the vertical interpolation
	662	mmean(i,j,l) = 0.
	663
	664	! 3.1 vertical interpolation
	665
	666	do n = 1,nspe
	667	call intrplf(pgcm,vmr,pinit,vmrinit(:,n),nalt)
	668	vmrint(n) = vmr
	669	! vmrint(n) = vmrinit(l,n)
	670	iq = niq(n)
	671	mmean(i,j,l) = mmean(i,j,l) + vmrint(n)*mmol(iq)
	672	end do
	673
	674	! 3.2 attribute mixing ratio: - all layers or only thermosphere
	675	! - with our without h2o
	676
	677
	678
	679	if (flagthermo == 0 .or. (flagthermo == 1 .and. exp(pgcm) < 0.1)) then
	680	do n = 1,nspe
	681	iq = niq(n)
	682	if (iq /= igcm_h2o_vap .or. flagh2o == 1) then
	683	pq(i,j,l,iq) = vmrint(n)*mmol(iq)/mmean(i,j,l)
	684	end if
	685	end do
	686	end if
	687
	688	end do
	689	end do
	690	end do
	691
	692	! set surface values of chemistry tracers to zero
	693
	694
	695	if (flagthermo == 0) then
	696	! NB: no problem for "surface water vapour" tracer which is always 0
	697	do n = 1,nspe
	698	iq = niq(n)
	699	qsurf(1:ngrid,iq) = 0.
	700	end do
	701	end if
	702
	703	! 3.3 initialization of tracers not contained in the initialization files
	704
	705	! methane : 10 ppbv
	706
	707	! if (igcm_ch4 /= 0) then
	708	! vmr = 10.e-9
	709	! do i = 1,nbp_lon+1
	710	! do j = 1,nbp_lat
	711	! do l = 1,nbp_lev
	712	! pq(i,j,l,igcm_ch4) = vmr*mmol(igcm_ch4)/mmean(i,j,l)
	713	! end do
	714	! end do
	715	! end do
	716	! ! set surface value to zero
	717	! qsurf(1:ngrid,igcm_ch4) = 0.
	718	! end if
	719
	720	! ions: 0
	721
	722
	723	! deallocations
	724
	725	deallocate(niq)
	726	deallocate(vmrinit)
	727	deallocate(vmrint)
	728
	729	end

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

Download in other formats:

Original Format