Context Navigation

hazecloud.F90 @ 3936

Last change on this file since 3936 was 3936, checked in by tbertrand, 2 months ago
PLUTO PCM : correcting a bug in hazecloud (wrong lyman alpha fluxes due to mu0 being negative during nighttime) + cleaning routines TB
File size: 12.6 KB

Line
1	subroutine hazecloud(ngrid,nlayer,nq,ptimestep,zday, &
2	pplay,pplev,zzlay,pq,pdq,pdist_sol,mu0,pfluxuv,zdqhaze, &
3	zdqphot_prec,zdqphot_ch4,zdqconv_prec,declin)
4	!zdqhaze_col)
5
6	use comgeomfi_h
7	use comcstfi_mod, only: pi, g
8	use tracer_h, only: igcm_haze, igcm_ch4_gas, igcm_n2, igcm_prec_haze, noms, mmol
9	use geometry_mod, only: longitude, latitude ! in radians
10	use callkeys_mod, only: hazeconservch4, haze_ch4proffix, diurnal, tcon_ch4, k_ch4, &
11	ncratio_ch4,triton,global1d
12	use datafile_mod, only: datadir
13	use write_output_mod, only: write_output
14
15	implicit none
16
17	!==================================================================
18	! Purpose
19	! -------
20	! Production of haze in the atmosphere
21	!
22	! Inputs
23	! ------
24	! ngrid Number of vertical columns
25	! nlayer Number of layers
26	! nq Number of tracers
27	! pplay(ngrid,nlayer) Pressure layers
28	! pplev(ngrid,nlayer+1) Pressure levels
29	! zzlay(ngrid,nlayer+1) mid layer altitude
30	! ptimestep Physical timestep
31	! zday, mu0 day clock, cos(incident flux angle)
32	! pq, pdq tracer MMR and tendency
33	! pdist_sol, declin Sun-planet distance, subsolar latitude
34	! pfluxuv Lyman alpha flux at Earth
35
36	! Outputs
37	! -------
38	! zdqhaze Tendency on tracers MMR haze
39	! zdqphot_prec Photolysis rate for haze precursors
40	! zdqphot_ch4 Photolysis rate for CH4
41	! zdqconv_prec Tendency for haze formation and precursor destruction
42
43	! Both
44	! ----
45	!
46	! Authors
47	! -------
48	! Tanguy Bertrand and Francois Forget (2014)
49	!
50	!==================================================================
51
52
53	!-----------------------------------------------------------------------
54	! Arguments
55
56	INTEGER ngrid, nlayer, nq
57	LOGICAL firstcall
58	SAVE firstcall
59	DATA firstcall/.true./
60	REAL ptimestep
61	REAL zday
62	REAL pplay(ngrid,nlayer),pplev(ngrid,nlayer+1)
63	REAL,INTENT(IN) :: zzlay(ngrid,nlayer) ! Mid-layer altitude
64	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)
65	REAL,INTENT(IN) :: pdq(ngrid,nlayer,nq)
66	REAL,INTENT(IN) :: pdist_sol ! distance SUN-pluto in AU
67	REAL,INTENT(IN) :: pfluxuv ! Lyman alpha flux at Earth at specific Ls (ph/cm2/s)
68	REAL,INTENT(IN) :: mu0(ngrid) ! cosinus of solar incident flux
69	REAL,INTENT(IN) :: declin ! Subsolar latitude
70	REAL,INTENT(OUT) :: zdqhaze(ngrid,nlayer,nq) ! Final tendancy
71	REAL,INTENT(OUT) :: zdqphot_ch4(ngrid,nlayer) ! tendancy due to photolysis of ch4
72	REAL,INTENT(OUT) :: zdqphot_prec(ngrid,nlayer) ! tendancy due to photolysis of ch4
73	REAL,INTENT(OUT) :: zdqconv_prec(ngrid,nlayer) ! tendancy due to conversion of
74	! prec_haze into haze
75	!-----------------------------------------------------------------------
76	! Local variables
77
78	INTEGER l,ig,iq
79	REAL zq_ch4(ngrid,nlayer)
80	REAL zq_prec(ngrid,nlayer)
81	REAL tauch4(nlayer)
82	REAL sigch4 ! aborption cross section ch4 cm-2 mol-1
83	REAL flym_sol_earth ! initial flux Earth ph.m-2.s-1
84	REAL flym_sol_pluto ! initial Lyman alpha SOLAR flux Pluto ph.m-2.s-1
85	REAL flym_ipm(ngrid) ! top of atm. Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1
86	REAL fluxlym_sol(nlayer+1) ! Local Solar flux Lyman alpha ph.m-2.s-1
87	REAL fluxlym_ipm(nlayer+1) ! Local IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1
88	REAL mu_ipm(ngrid) ! local Mean incident flux for IPM Lyman alpha photons
89	REAL mu_sol(ngrid) ! local Mean incident flux for solar lyman alpha photons
90
91	REAL gradflux(nlayer) ! gradient flux Lyman alpha ph.m-2.s-1
92	REAL kch4 ! fraction of Carbon from ch4 directly dissociated
93	! by prec_haze
94	CHARACTER(len=10) :: tname
95	REAL ncratio ! ration Nitrogen/Carbon observed in tholins
96	REAL tcon ! Time constant: conversion in aerosol
97	REAL avogadro ! avogadro constant
98
99	integer Nfine,ifine
100	parameter(Nfine=600)
101	character(len=100) :: file_path
102	real,save :: levdat(Nfine),vmrdat(Nfine)
103	real :: vmrch4(ngrid,nlayer) ! vmr ch4 from vmrch4_proffix
104
105	REAL longit
106	REAL latit
107	REAL valmin
108	REAL valmax
109	REAL valmin_dl
110	REAL puis
111	REAL dist
112	REAL long2
113
114	! Diagnostics (temporary)
115	REAL fluxlym_sol_tot(ngrid)
116	REAL fluxlym_ipm_tot(ngrid)
117
118	!-----------------------------------------------------------------------
119	! Input parameters
120
121	avogadro = 6.022141e23
122	sigch4 = 1.85e-17 ! aborption cross section of ch4 in cm-2 mol-1
123	! Initial Solar flux Lyman alpha at Earth (1AU)
124	flym_sol_earth=pfluxuv*1.e15 ! ph.m-2.s-1 -> about 5e+11 ph.cm-2.s-1
125	! Initial Solar flux Lyman alpha on Pluto
126	flym_sol_pluto=flym_sol_earth/(pdist_sol*pdist_sol) ! ph.m-2.s-1
127	! option decrease by 12% the initial Solar Lyman alpha flux to account for Interplanetary H absorption:
128	! Cf Fig 3 Gladstone et al 2014 : Lyalpha at Pluto
129	flym_sol_pluto=flym_sol_pluto*0.878
130
131	! Time constant of conversion in aerosol [second]
132	! To be explored: 1.E5 - 1.E9
133	tcon= tcon_ch4
134	! Parameter of conversion precurseur to haze
135	kch4=k_ch4
136	ncratio=ncratio_ch4 ! boost for haze considering nitrogen contribution
137	! ratio n/c : =0.25 if there is 1N for 3C
138	IF (firstcall) then
139	write(,) 'hazecloud: haze parameters:'
140	write(,) 'tcon, kch4, ncratio = ' , tcon, kch4, ncratio
141	firstcall=.false.
142	ENDIF
143	! note: mu0 = cos(solar zenith angle)
144	! max(mu0) at lat = declin
145
146	!-----------------------------------------------------------------------
147	! Total IPM Lyman alpha flux at top of atmosphere
148
149	! Top of atm. Incident IPM (Interplanetary) flux Lyman alpha ph.m-2.s-1
150	! fit Fig. 4 in Randall Gladstone et al. (2015)
151	! Integration over semi sphere of Gladstone data.
152	! 145 R averaged over the sky
153	! 72.5 R in average over the visible hemisphere
154	IF (ngrid.eq.1.and..not.global1d) THEN
155	mu_sol(1) = mu0(1)
156	mu_ipm(1) = max(mu_sol(1), 0.5)
157	flym_ipm(1)=mu0(1)*72.5e10
158	ELSE IF (ngrid.eq.1.and.global1d) THEN
159	mu_sol(1) = mu0(1) ! Full visible disk
160	mu_ipm(1) = mu0(1)/2. ! Full sphere (day+night)
161	flym_ipm(1)=mu_ipm(1)*145.e10
162	! Proportional to lyman alpha solar flux (reference 2015 : 5e11 ph/cm2/s) ?
163	! In theory, IPM flux remains relatively constant further than 15 AU
164	flym_ipm(1)=flym_ipm(1)*pfluxuv/5.
165	ELSE IF (.not.diurnal) THEN
166	mu_sol(:) = max(mu0(:),0.)
167	mu_ipm(:) = max(mu_sol(:), 0.5)
168	flym_ipm(:)= mu_sol(:)*72.5e10
169	ELSE ! case with full fit to Gladstone et al. results
170	! 1) get longitude/latitude (in radian) of anti-subsolar point (max de mu0 - 180)
171	longit=-(zday-int(zday))2.pi
172	IF (longit.LE.-pi) THEN
173	longit=longit+2.*pi
174	ENDIF
175	latit=-declin ! rad
176	! 2) Define input parameter for the fit
177	valmin=48.74e10 ! minimum value of ipm flux in ph/m2/s
178	valmax=115.014e10 ! maximal value of ipm flux in ph/m2/s
179	valmin_dl=74.5e10 ! daylight minimum value
180	puis=3.5
181	! 3) Loop for each location and fit
182	DO ig=1,ngrid
183	! calculation of cosinus of incident angle for IPM flux
184	mu_sol(ig) = max(mu0(ig),0.)
185	mu_ipm(ig) = max(mu_sol(ig), 0.5)
186	IF (mu_sol(ig).LE.0.) THEN ! Nighttime
187	! Distance to subsolar point
188	dist=acos(sin(latitude(ig))sin(latit)+cos(latitude(ig)) &
189	cos(latit)cos(longit-longitude(ig)))180/pi
190	IF (dist.gt.90) dist=90
191	flym_ipm(ig)=(valmin_dl-valmin)/(90.*puis)(dist)**puis &
192	+valmin
193	ELSE ! Daytime
194	flym_ipm(ig)= mu_sol(ig)*(valmax-valmin_dl)+valmin_dl
195	ENDIF
196	! Proportional to lyman alpha solar flux (reference 2015 : 5e11 ph/cm2/s) ?
197	! In theory, IPM flux remains relatively constant further than 15 AU
198	! flym_ipm(ig)=flym_ipm(ig)*pfluxuv/5.
199	ENDDO ! nlayer
200
201	ENDIF ! ngrid
202
203	!-----------------------------------------------------------------------
204	! Methane profile
205
206	! If fixed profile of CH4 gas
207	if (haze_ch4proffix) then
208	if (triton) then
209	file_path=trim(datadir)//'/gas_prop/vmr_ch4_triton.txt'
210	else
211	file_path=trim(datadir)//'/gas_prop/vmr_ch4.txt'
212	endif
213	open(115,file=file_path,form='formatted')
214	do ifine=1,Nfine
215	read(115,*) levdat(ifine), vmrdat(ifine)
216	enddo
217	close(115)
218	endif
219	! Initialization:
220	vmrch4(:,:)=0.
221
222	! Diagnostics (temporary)
223	fluxlym_sol_tot(:)=flym_sol_pluto*mu_sol(:)
224	fluxlym_ipm_tot(:)=flym_ipm(:)
225	call write_output("fluxlym_sol","solar lyman alpha flux","",fluxlym_sol_tot)
226	call write_output("fluxlym_ipm","IPM lyman alpha flux","",fluxlym_ipm_tot)
227	!-----------------------------------------------------------------------
228	! Main Loop
229
230	DO ig=1,ngrid
231
232	!! Initializations
233	zq_ch4(ig,:)=0.
234	zq_prec(ig,:)=0.
235	zdqconv_prec(ig,:)=0.
236	zdqhaze(ig,:,:)=0.
237	zdqphot_prec(ig,:)=0.
238	zdqphot_ch4(ig,:)=0.
239	tauch4(:)=0.
240	gradflux(:)=0.
241	fluxlym_sol(1:nlayer)=0.
242
243	fluxlym_sol(nlayer+1)=flym_sol_pluto*mu_sol(ig) !
244	fluxlym_ipm(nlayer+1)= flym_ipm(ig)
245
246	!! Interpolate CH4 MMR on the model vertical grid
247	if (haze_ch4proffix) then
248	CALL interp_line(levdat,vmrdat,Nfine, &
249	zzlay(ig,:)/1000.,vmrch4(ig,:),nlayer)
250	endif
251
252	DO l=nlayer,1,-1
253
254	!! Actualisation tracer ch4 and prec_haze
255	if (haze_ch4proffix) then
256	zq_ch4(ig,l)=vmrch4(ig,l)*mmol(igcm_ch4_gas) / &
257	(100.*mmol(igcm_n2))
258	else
259	zq_ch4(ig,l)=pq(ig,l,igcm_ch4_gas) + &
260	pdq(ig,l,igcm_ch4_gas)*ptimestep
261	endif
262
263	zq_prec(ig,l)=pq(ig,l,igcm_prec_haze)+ &
264	pdq(ig,l,igcm_prec_haze)*ptimestep
265
266	!! Sanity check
267	if (zq_ch4(ig,l).lt.0.) then
268	zq_ch4(ig,l)=0.
269	endif
270	if (zq_prec(ig,l).lt.0.) then
271	zq_prec(ig,l)=0.
272	endif
273
274	!! Calculation of CH4 optical depth in Lyman alpha for each layer l
275	tauch4(l)=sigch41.e-4avogadro* &
276	(zq_ch4(ig,l)/(mmol(igcm_ch4_gas)1.e-3)) &
277	(pplev(ig,l)-pplev(ig,l+1))/g
278
279	!! Calculation of Flux in each layer l
280	if (mu_sol(ig).gt.1.e-6) then
281	fluxlym_sol(l)=fluxlym_sol(l+1)*exp(-tauch4(l)/mu_sol(ig))
282	endif
283
284	fluxlym_ipm(l)=fluxlym_ipm(l+1)*exp(-tauch4(l)/mu_ipm(ig))
285
286	gradflux(l)=fluxlym_sol(l+1)-fluxlym_sol(l) &
287	+ fluxlym_ipm(l+1)-fluxlym_ipm(l)
288
289	!! tendancy on ch4
290	!! considering 1 photon destroys 1 ch4 by photolysis
291	zdqphot_ch4(ig,l)=-mmol(igcm_ch4_gas)1.e-3/avogadro &
292	gradflux(l)*g/(pplev(ig,l)-pplev(ig,l+1))
293
294	!! tendency of precursors created by photolysis of ch4
295	zdqphot_prec(ig,l)=-zdqphot_ch4(ig,l)
296
297	!! update of precursors mass: zq_prec
298	zq_prec(ig,l)=zq_prec(ig,l)+ &
299	zdqphot_prec(ig,l)*ptimestep
300
301	!! Conversion of prec_haze into haze
302	!! controlled by the time constant tcon
303	!! New tendancy for prec_haze
304	zdqconv_prec(ig,l) = -zq_prec(ig,l)* &
305	(1-exp(-ptimestep/tcon))/ptimestep
306
307	ENDDO ! nlayer
308
309	!! Final tendancy for prec_haze, haze and CH4
310	DO iq=1,nq
311	tname=noms(iq)
312	if (tname(1:4).eq."haze") then
313	zdqhaze(ig,:,iq) = -zdqconv_prec(ig,:)* &
314	kch4(1.+ncratio14./12.)
315	else if (noms(iq).eq."prec_haze") then
316	zdqhaze(ig,:,igcm_prec_haze)= zdqphot_prec(ig,:) + &
317	zdqconv_prec(ig,:)
318	else if (noms(iq).eq."ch4_gas".and.(.not.hazeconservch4).and.(.not.haze_ch4proffix)) then
319	zdqhaze(ig,:,igcm_ch4_gas)= zdqphot_ch4(ig,:)
320	endif
321	ENDDO
322
323	ENDDO ! ngrid
324
325	end subroutine hazecloud
326

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

Download in other formats:

Original Format