Context Navigation

inichim_newstart.F90 @ 3807

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

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

Download in other formats:

Original Format