Context Navigation

inichim_newstart.F90 @ 3567

Last change on this file since 3567 was 3466, checked in by emillour, 2 months ago

Mars PCM:
More tidying in aeronomars:

remove unused "inv.F" and remove "dtridgl.F" which is not used here and is a duplicate of the "dtridgl" routine in phymars/swr_toon.F
turn aeronomars routines to modules, for those which aren't in modules yet.

File size: 24.7 KB

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

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

Download in other formats:

Original Format