Context Navigation

initracer.F @ 1892

Last change on this file since 1892 was 1801, checked in by bclmd, 8 years ago
Adding photochemistry to LMDZ Generic
File size: 16.4 KB

Rev	Line
[787]	1	SUBROUTINE initracer(ngrid,nq,nametrac)
[135]	2
[787]	3	use surfdat_h
	4	USE tracer_h
[1397]	5	USE callkeys_mod, only: water
[135]	6	IMPLICIT NONE
	7	c=======================================================================
	8	c subject:
	9	c --------
	10	c Initialization related to tracer
	11	c (transported dust, water, chemical species, ice...)
	12	c
	13	c Name of the tracer
	14	c
	15	c Test of dimension :
	16	c Initialize COMMON tracer in tracer.h, using tracer names provided
[787]	17	c by the argument nametrac
[135]	18	c
	19	c author: F.Forget
	20	c ------
	21	c Ehouarn Millour (oct. 2008) identify tracers by their names
	22	c=======================================================================
	23
[787]	24	integer :: ngrid,nq
	25
	26	! real qsurf(ngrid,nq) ! tracer on surface (e.g. kg.m-2)
	27	! real co2ice(ngrid) ! co2 ice mass on surface (e.g. kg.m-2)
[135]	28	character(len=20) :: txt ! to store some text
	29	integer iq,ig,count
	30	real r0_lift , reff_lift
	31	! logical :: oldnames ! =.true. if old tracer naming convention (q01,...)
	32
[787]	33	character*20 nametrac(nq) ! name of the tracer from dynamics
[135]	34
[787]	35
[135]	36	c-----------------------------------------------------------------------
[787]	37	c radius(nq) ! aerosol particle radius (m)
	38	c rho_q(nq) ! tracer densities (kg.m-3)
	39	c qext(nq) ! Single Scat. Extinction coeff at 0.67 um
	40	c alpha_lift(nq) ! saltation vertical flux/horiz flux ratio (m-1)
	41	c alpha_devil(nq) ! lifting coeeficient by dust devil
[135]	42	c rho_dust ! Mars dust density
	43	c rho_ice ! Water ice density
	44	c doubleq ! if method with mass (iq=1) and number(iq=2) mixing ratio
	45	c varian ! Characteristic variance of log-normal distribution
	46	c-----------------------------------------------------------------------
	47
[1621]	48	nqtot=nq
[787]	49	!! we allocate once for all arrays in common in tracer_h.F90
	50	!! (supposedly those are not used before call to initracer)
[1801]	51	IF (.NOT.ALLOCATED(noms)) ALLOCATE(noms(nq))
	52	IF (.NOT.ALLOCATED(mmol)) ALLOCATE(mmol(nq))
	53	IF (.NOT.ALLOCATED(radius)) ALLOCATE(radius(nq))
	54	IF (.NOT.ALLOCATED(rho_q)) ALLOCATE(rho_q(nq))
	55	IF (.NOT.ALLOCATED(qext)) ALLOCATE(qext(nq))
	56	IF (.NOT.ALLOCATED(alpha_lift)) ALLOCATE(alpha_lift(nq))
	57	IF (.NOT.ALLOCATED(alpha_devil)) ALLOCATE(alpha_devil(nq))
	58	IF (.NOT.ALLOCATED(qextrhor)) ALLOCATE(qextrhor(nq))
	59	IF (.NOT.ALLOCATED(igcm_dustbin)) ALLOCATE(igcm_dustbin(nq))
[787]	60	!! initialization
	61	alpha_lift(:)=0.
	62	alpha_devil(:)=0.
[1764]	63
	64	! Added by JVO 2017 : these arrays are handled later
	65	! -> initialization is the least we can do, please !!!
	66	radius(:)=0.
	67	qext(:)=0.
[787]	68
[135]	69	! Initialization: get tracer names from the dynamics and check if we are
	70	! using 'old' tracer convention ('q01',q02',...)
	71	! or new convention (full tracer names)
	72	! check if tracers have 'old' names
	73
	74	! copy tracer names from dynamics
[787]	75	do iq=1,nq
	76	noms(iq)=nametrac(iq)
[135]	77	enddo
	78
	79
	80	! Identify tracers by their names: (and set corresponding values of mmol)
	81	! 0. initialize tracer indexes to zero:
	82	! NB: igcm_* indexes are commons in 'tracer.h'
[787]	83	do iq=1,nq
[135]	84	igcm_dustbin(iq)=0
	85	enddo
	86	igcm_dust_mass=0
	87	igcm_dust_number=0
	88	igcm_h2o_vap=0
	89	igcm_h2o_ice=0
	90	igcm_co2=0
	91	igcm_co=0
	92	igcm_o=0
	93	igcm_o1d=0
	94	igcm_o2=0
	95	igcm_o3=0
	96	igcm_h=0
	97	igcm_h2=0
	98	igcm_oh=0
	99	igcm_ho2=0
	100	igcm_h2o2=0
	101	igcm_n2=0
[1801]	102	igcm_n=0
	103	igcm_n2d=0
	104	igcm_no=0
	105	igcm_no2=0
[135]	106	igcm_ar=0
	107	igcm_ar_n2=0
	108	igcm_co2_ice=0
	109
[1801]	110	igcm_ch4=0
	111	igcm_ch3=0
	112	igcm_ch=0
	113	igcm_3ch2=0
	114	igcm_1ch2=0
	115	igcm_cho=0
	116	igcm_ch2o=0
	117	igcm_ch3o=0
	118	igcm_c=0
	119	igcm_c2=0
	120	igcm_c2h=0
	121	igcm_c2h2=0
	122	igcm_c2h3=0
	123	igcm_c2h4=0
	124	igcm_c2h6=0
	125	igcm_ch2co=0
	126	igcm_ch3co=0
	127	igcm_hcaer=0
	128
	129
[135]	130	write(,) 'initracer: noms() ', noms
	131
	132
[253]	133	!print*,'Setting dustbin = 0 in initracer.F'
	134	!dustbin=0
[135]	135
	136	! 1. find dust tracers
	137	count=0
	138	! if (dustbin.gt.0) then
[787]	139	! do iq=1,nq
[135]	140	! txt=" "
	141	! write(txt,'(a4,i2.2)')'dust',count+1
	142	! if (noms(iq).eq.txt) then
	143	! count=count+1
	144	! igcm_dustbin(count)=iq
	145	! mmol(iq)=100.
	146	! endif
[787]	147	! enddo !do iq=1,nq
[135]	148	! endif ! of if (dustbin.gt.0)
	149
	150
	151	! if (doubleq) then
[787]	152	! do iq=1,nq
[135]	153	! if (noms(iq).eq."dust_mass") then
	154	! igcm_dust_mass=iq
	155	! count=count+1
	156	! endif
	157	! if (noms(iq).eq."dust_number") then
	158	! igcm_dust_number=iq
	159	! count=count+1
	160	! endif
	161	! enddo
	162	! endif ! of if (doubleq)
	163	! 2. find chemistry and water tracers
[787]	164	do iq=1,nq
[135]	165	if (noms(iq).eq."co2") then
	166	igcm_co2=iq
	167	mmol(igcm_co2)=44.
	168	count=count+1
	169	! write(,) 'co2: count=',count
	170	endif
	171	if (noms(iq).eq."co2_ice") then
	172	igcm_co2_ice=iq
	173	mmol(igcm_co2_ice)=44.
	174	count=count+1
	175	! write(,) 'co2_ice: count=',count
	176	endif
	177	if (noms(iq).eq."h2o_vap") then
	178	igcm_h2o_vap=iq
	179	mmol(igcm_h2o_vap)=18.
	180	count=count+1
	181	! write(,) 'h2o_vap: count=',count
	182	endif
	183	if (noms(iq).eq."h2o_ice") then
	184	igcm_h2o_ice=iq
	185	mmol(igcm_h2o_ice)=18.
	186	count=count+1
	187	! write(,) 'h2o_ice: count=',count
	188	endif
[1801]	189	if (noms(iq).eq."co") then
	190	igcm_co=iq
	191	mmol(igcm_co)=28.
	192	count=count+1
	193	endif
	194	if (noms(iq).eq."o") then
	195	igcm_o=iq
	196	mmol(igcm_o)=16.
	197	count=count+1
	198	endif
	199	if (noms(iq).eq."o1d") then
	200	igcm_o1d=iq
	201	mmol(igcm_o1d)=16.
	202	count=count+1
	203	endif
	204	if (noms(iq).eq."o2") then
	205	igcm_o2=iq
	206	mmol(igcm_o2)=32.
	207	count=count+1
	208	endif
	209	if (noms(iq).eq."o3") then
	210	igcm_o3=iq
	211	mmol(igcm_o3)=48.
	212	count=count+1
	213	endif
	214	if (noms(iq).eq."h") then
	215	igcm_h=iq
	216	mmol(igcm_h)=1.
	217	count=count+1
	218	endif
	219	if (noms(iq).eq."h2") then
	220	igcm_h2=iq
	221	mmol(igcm_h2)=2.
	222	count=count+1
	223	endif
	224	if (noms(iq).eq."oh") then
	225	igcm_oh=iq
	226	mmol(igcm_oh)=17.
	227	count=count+1
	228	endif
	229	if (noms(iq).eq."ho2") then
	230	igcm_ho2=iq
	231	mmol(igcm_ho2)=33.
	232	count=count+1
	233	endif
	234	if (noms(iq).eq."h2o2") then
	235	igcm_h2o2=iq
	236	mmol(igcm_h2o2)=34.
	237	count=count+1
	238	endif
	239	if (noms(iq).eq."n2") then
	240	igcm_n2=iq
	241	mmol(igcm_n2)=28.
	242	count=count+1
	243	endif
	244	if (noms(iq).eq."ch4") then
	245	igcm_ch4=iq
	246	mmol(igcm_ch4)=16.
	247	count=count+1
	248	endif
	249	if (noms(iq).eq."ar") then
	250	igcm_ar=iq
	251	mmol(igcm_ar)=40.
	252	count=count+1
	253	endif
	254	if (noms(iq).eq."n") then
	255	igcm_n=iq
	256	mmol(igcm_n)=14.
	257	count=count+1
	258	endif
	259	if (noms(iq).eq."no") then
	260	igcm_no=iq
	261	mmol(igcm_no)=30.
	262	count=count+1
	263	endif
	264	if (noms(iq).eq."no2") then
	265	igcm_no2=iq
	266	mmol(igcm_no2)=46.
	267	count=count+1
	268	endif
	269	if (noms(iq).eq."n2d") then
	270	igcm_n2d=iq
	271	mmol(igcm_n2d)=28.
	272	count=count+1
	273	endif
	274	if (noms(iq).eq."ch3") then
	275	igcm_ch3=iq
	276	mmol(igcm_ch3)=15.
	277	count=count+1
	278	endif
	279	if (noms(iq).eq."ch") then
	280	igcm_ch=iq
	281	mmol(igcm_ch)=13.
	282	count=count+1
	283	endif
	284	if (noms(iq).eq."3ch2") then
	285	igcm_3ch2=iq
	286	mmol(igcm_3ch2)=14.
	287	count=count+1
	288	endif
	289	if (noms(iq).eq."1ch2") then
	290	igcm_1ch2=iq
	291	mmol(igcm_1ch2)=14.
	292	count=count+1
	293	endif
	294	if (noms(iq).eq."cho") then
	295	igcm_cho=iq
	296	mmol(igcm_cho)=29.
	297	count=count+1
	298	endif
	299	if (noms(iq).eq."ch2o") then
	300	igcm_ch2o=iq
	301	mmol(igcm_ch2o)=30.
	302	count=count+1
	303	endif
	304	if (noms(iq).eq."ch3o") then
	305	igcm_ch3o=iq
	306	mmol(igcm_ch3o)=31.
	307	count=count+1
	308	endif
	309	if (noms(iq).eq."c") then
	310	igcm_c=iq
	311	mmol(igcm_c)=12.
	312	count=count+1
	313	endif
	314	if (noms(iq).eq."c2") then
	315	igcm_c2=iq
	316	mmol(igcm_c2)=24.
	317	count=count+1
	318	endif
	319	if (noms(iq).eq."c2h") then
	320	igcm_c2h=iq
	321	mmol(igcm_c2h)=25.
	322	count=count+1
	323	endif
	324	if (noms(iq).eq."c2h2") then
	325	igcm_c2h2=iq
	326	mmol(igcm_c2h2)=26.
	327	count=count+1
	328	endif
	329	if (noms(iq).eq."c2h3") then
	330	igcm_c2h3=iq
	331	mmol(igcm_c2h3)=27.
	332	count=count+1
	333	endif
	334	if (noms(iq).eq."c2h4") then
	335	igcm_c2h4=iq
	336	mmol(igcm_c2h4)=28.
	337	count=count+1
	338	endif
	339	if (noms(iq).eq."c2h6") then
	340	igcm_c2h6=iq
	341	mmol(igcm_c2h6)=30.
	342	count=count+1
	343	endif
	344	if (noms(iq).eq."ch2co") then
	345	igcm_ch2co=iq
	346	mmol(igcm_ch2co)=42.
	347	count=count+1
	348	endif
	349	if (noms(iq).eq."ch3co") then
	350	igcm_ch3co=iq
	351	mmol(igcm_ch3co)=43.
	352	count=count+1
	353	endif
	354	if (noms(iq).eq."hcaer") then
	355	igcm_hcaer=iq
	356	mmol(igcm_hcaer)=50.
	357	count=count+1
	358	endif
[787]	359	enddo ! of do iq=1,nq
[135]	360
	361	! check that we identified all tracers:
[787]	362	if (count.ne.nq) then
[135]	363	write(,) "initracer: found only ",count," tracers"
[787]	364	write(,) " expected ",nq
[135]	365	do iq=1,count
	366	write(,)' ',iq,' ',trim(noms(iq))
	367	enddo
[1297]	368	! stop
[135]	369	else
	370	write(,) "initracer: found all expected tracers, namely:"
[787]	371	do iq=1,nq
[135]	372	write(,)' ',iq,' ',trim(noms(iq))
	373	enddo
	374	endif
	375
	376
	377	c------------------------------------------------------------
	378	c Initialisation tracers ....
	379	c------------------------------------------------------------
[787]	380	call zerophys(nq,rho_q)
[135]	381
	382	rho_dust=2500. ! Mars dust density (kg.m-3)
	383	rho_ice=920. ! Water ice density (kg.m-3)
	384	rho_co2=1620. ! CO2 ice density (kg.m-3)
	385
	386
	387
	388	c$$$ if (doubleq) then
	389	c$$$c "doubleq" technique
	390	c$$$c -------------------
	391	c$$$c (transport of mass and number mixing ratio)
	392	c$$$c iq=1: Q mass mixing ratio, iq=2: N number mixing ratio
	393	c$$$
[787]	394	c$$$ if( (nq.lt.2).or.(water.and.(nq.lt.3)) ) then
	395	c$$$ write(,)'initracer: nq is too low : nq=', nq
[135]	396	c$$$ write(,)'water= ',water,' doubleq= ',doubleq
	397	c$$$ end if
	398	c$$$
	399	c$$$ varian=0.637 ! Characteristic variance
	400	c$$$ qext(igcm_dust_mass)=3.04 ! reference extinction at 0.67 um for ref dust
	401	c$$$ qext(igcm_dust_number)=3.04 ! reference extinction at 0.67 um for ref dust
	402	c$$$ rho_q(igcm_dust_mass)=rho_dust
	403	c$$$ rho_q(igcm_dust_number)=rho_dust
	404	c$$$
	405	c$$$c Intermediate calcul for computing geometric mean radius r0
	406	c$$$c as a function of mass and number mixing ratio Q and N
	407	c$$$c (r0 = (r3n_q * Q/ N)
	408	c$$$ r3n_q = exp(-4.5varian2)(3./4.)/(pi*rho_dust)
	409	c$$$
	410	c$$$c Intermediate calcul for computing effective radius reff
	411	c$$$c from geometric mean radius r0
	412	c$$$c (reff = ref_r0 * r0)
	413	c$$$ ref_r0 = exp(2.5varian*2)
	414	c$$$
	415	c$$$c lifted dust :
	416	c$$$c '''''''''''
	417	c$$$ reff_lift = 3.e-6 ! Effective radius of lifted dust (m)
	418	c$$$ alpha_devil(igcm_dust_mass)=9.e-9 ! dust devil lift mass coeff
	419	c$$$ alpha_lift(igcm_dust_mass)=3.0e-15 ! Lifted mass coeff
	420	c$$$
	421	c$$$ r0_lift = reff_lift/ref_r0
	422	c$$$ alpha_devil(igcm_dust_number)=r3n_q*
	423	c$$$ & alpha_devil(igcm_dust_mass)/r0_lift**3
	424	c$$$ alpha_lift(igcm_dust_number)=r3n_q*
	425	c$$$ & alpha_lift(igcm_dust_mass)/r0_lift**3
	426	c$$$
	427	c$$$c Not used:
	428	c$$$ radius(igcm_dust_mass) = 0.
	429	c$$$ radius(igcm_dust_number) = 0.
	430	c$$$
	431	c$$$ else
	432
	433
[253]	434	c$$$ if (dustbin.gt.1) then
	435	c$$$ print*,'ATTENTION:',
	436	c$$$ $ ' properties of dust need input in initracer !!!'
	437	c$$$ stop
	438	c$$$
	439	c$$$ else if (dustbin.eq.1) then
	440	c$$$
	441	c$$$c This will be used for 1 dust particle size:
	442	c$$$c ------------------------------------------
	443	c$$$ radius(igcm_dustbin(1))=3.e-6
	444	c$$$ Qext(igcm_dustbin(1))=3.04
	445	c$$$ alpha_lift(igcm_dustbin(1))=0.0e-6
	446	c$$$ alpha_devil(igcm_dustbin(1))=7.65e-9
	447	c$$$ qextrhor(igcm_dustbin(1))=(3./4.)*Qext(igcm_dustbin(1))
	448	c$$$ & /(rho_dust*radius(igcm_dustbin(1)))
	449	c$$$ rho_q(igcm_dustbin(1))=rho_dust
	450	c$$$
	451	c$$$ endif
	452	c$$$! end if ! (doubleq)
[135]	453
	454	c Initialization for water vapor
	455	c ------------------------------
	456	if(water) then
	457	radius(igcm_h2o_vap)=0.
	458	Qext(igcm_h2o_vap)=0.
	459	alpha_lift(igcm_h2o_vap) =0.
	460	alpha_devil(igcm_h2o_vap)=0.
	461	qextrhor(igcm_h2o_vap)= 0.
	462
	463	c "Dryness coefficient" controlling the evaporation and
	464	c sublimation from the ground water ice (close to 1)
	465	c HERE, the goal is to correct for the fact
	466	c that the simulated permanent water ice polar caps
	467	c is larger than the actual cap and the atmospheric
	468	c opacity not always realistic.
	469
[253]	470
[787]	471	! if(ngrid.eq.1)
[253]	472
	473
	474	! to be modified for BC+ version?
[787]	475
	476	!! this is defined in surfdat_h.F90
	477	IF (.not.ALLOCATED(dryness)) ALLOCATE(dryness(ngrid))
	478	IF (.not.ALLOCATED(watercaptag)) ALLOCATE(watercaptag(ngrid))
	479
	480	do ig=1,ngrid
	481	if (ngrid.ne.1) watercaptag(ig)=.false.
[135]	482	dryness(ig) = 1.
	483	enddo
	484
[253]	485
	486
	487
[135]	488	! IF (caps) THEN
	489	c Perennial H20 north cap defined by watercaptag=true (allows surface to be
	490	c hollowed by sublimation in vdifc).
[787]	491	! do ig=1,ngrid
[135]	492	! if (lati(ig)*180./pi.gt.84) then
[787]	493	! if (ngrid.ne.1) watercaptag(ig)=.true.
[135]	494	! dryness(ig) = 1.
	495	c Use the following cap definition for high spatial resolution (latitudinal bin <= 5 deg)
	496	c if (lati(ig)*180./pi.lt.85.and.long(ig).ge.0) then
[787]	497	c if (ngrid.ne.1) watercaptag(ig)=.true.
[135]	498	c dryness(ig) = 1.
	499	c endif
	500	c if (lati(ig)*180./pi.ge.85) then
[787]	501	c if (ngrid.ne.1) watercaptag(ig)=.true.
[135]	502	c dryness(ig) = 1.
	503	c endif
	504	! endif ! (lati>80 deg)
[787]	505	! end do ! (ngrid)
[135]	506	! ENDIF ! (caps)
	507
[787]	508	! if(iceparty.and.(nq.ge.2)) then
[135]	509
	510	radius(igcm_h2o_ice)=3.e-6
	511	rho_q(igcm_h2o_ice)=rho_ice
	512	Qext(igcm_h2o_ice)=0.
	513	! alpha_lift(igcm_h2o_ice) =0.
	514	! alpha_devil(igcm_h2o_ice)=0.
	515	qextrhor(igcm_h2o_ice)= (3./4.)*Qext(igcm_h2o_ice)
	516	$ / (rho_ice*radius(igcm_h2o_ice))
	517
	518
	519
[787]	520	! elseif(iceparty.and.(nq.lt.2)) then
	521	! write(,) 'nq is too low : nq=', nq
[135]	522	! write(,) 'water= ',water,' iceparty= ',iceparty
	523	! endif
	524
	525	end if ! (water)
	526
[1801]	527
	528	!
	529	! some extra (possibly redundant) sanity checks for tracers:
	530	! ---------------------------------------------------------
	531	if (water) then
	532	! verify that we indeed have h2o_vap and h2o_ice tracers
	533	if (igcm_h2o_vap.eq.0) then
	534	write(,) "initracer: error !!"
	535	write(,) " cannot use water option without ",
	536	& "an h2o_vap tracer !"
	537	stop
	538	endif
	539	if (igcm_h2o_ice.eq.0) then
	540	write(,) "initracer: error !!"
	541	write(,) " cannot use water option without ",
	542	& "an h2o_ice tracer !"
	543	stop
	544	endif
	545	endif
	546
	547
[135]	548	c Output for records:
	549	c ~~~~~~~~~~~~~~~~~~
	550	write(,)
	551	Write(,) '******** initracer : dust transport parameters :'
	552	write(,) 'alpha_lift = ', alpha_lift
	553	write(,) 'alpha_devil = ', alpha_devil
	554	write(,) 'radius = ', radius
	555	! if(doubleq) then
	556	! write(,) 'reff_lift (um) = ', reff_lift
	557	! write(,) 'size distribution variance = ', varian
	558	! write(,) 'r3n_q , ref_r0 : ', r3n_q , ref_r0
	559	! end if
	560	write(,) 'Qext = ', qext
	561	write(,)
	562
	563	end

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

Download in other formats:

Original Format