Context Navigation

inichim_newstart.F90 @ 2613

Last change on this file since 2613 was 2579, checked in by lrossi, 3 years ago

29/10/2021 == LR

MARS GCM
Fixing missing case for hdo_ice in inichim_newstart
Should fix issue #25 on gitlab

File size: 22.3 KB

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

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

Download in other formats:

Original Format