Context Navigation

inichim_newstart.F90 @ 3026

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

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

Download in other formats:

Original Format