Context Navigation

initracer.F @ 2176

Last change on this file since 2176 was 1980, checked in by emillour, 6 years ago
Generic physics: Follow-up of change in maximum length of tracer names in the dynamics EM
File size: 16.2 KB

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

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

Download in other formats:

Original Format