Context Navigation

aeropacity.F90 @ 1974

Last change on this file since 1974 was 1677, checked in by sglmd, 8 years ago

Two aerosol kinds added for giant planets: one is a compact (NH3) cloud where the opacity, particle size and bottom pressure level are taken as inputs (a scale height of 0.2 is hard-coded to simulate a compact cloud). Corresponds to option aeronh3. The other one is not generic at all and corresponds to the auroral, stratospheric aerosols observed on Jupiter...(option aeroaurora=.false. by default)

File size: 20.1 KB

Line
1	Subroutine aeropacity(ngrid,nlayer,nq,pplay,pplev,pq, &
2	aerosol,reffrad,QREFvis3d,QREFir3d,tau_col,cloudfrac,totcloudfrac,clearsky)
3
4	use radinc_h, only : L_TAUMAX,naerkind
5	use aerosol_mod
6	USE tracer_h, only: noms,rho_co2,rho_ice
7	use comcstfi_mod, only: g, pi
8	use geometry_mod, only: latitude
9	use callkeys_mod, only: aerofixco2,aerofixh2o,kastprof,cloudlvl, &
10	CLFvarying,CLFfixval,dusttau, &
11	pres_bottom_tropo,pres_top_tropo,obs_tau_col_tropo, &
12	pres_bottom_strato,pres_top_strato,obs_tau_col_strato, &
13	tau_nh3_cloud, pres_nh3_cloud
14
15	implicit none
16
17	!==================================================================
18	!
19	! Purpose
20	! -------
21	! Compute aerosol optical depth in each gridbox.
22	!
23	! Authors
24	! -------
25	! F. Forget
26	! F. Montmessin (water ice scheme)
27	! update J.-B. Madeleine (2008)
28	! dust removal, simplification by Robin Wordsworth (2009)
29	!
30	! Input
31	! -----
32	! ngrid Number of horizontal gridpoints
33	! nlayer Number of layers
34	! nq Number of tracers
35	! pplev Pressure (Pa) at each layer boundary
36	! pq Aerosol mixing ratio
37	! reffrad(ngrid,nlayer,naerkind) Aerosol effective radius
38	! QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients
39	! QREFir3d(ngrid,nlayer,naerkind) / at reference wavelengths
40	!
41	! Output
42	! ------
43	! aerosol Aerosol optical depth in layer l, grid point ig
44	! tau_col Total column optical depth at grid point ig
45	!
46	!=======================================================================
47
48	INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
49	INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
50	INTEGER,INTENT(IN) :: nq ! number of tracers
51	REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
52	REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
53	REAL,INTENT(IN) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
54	REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! aerosol optical depth
55	REAL,INTENT(IN) :: reffrad(ngrid,nlayer,naerkind) ! aerosol effective radius
56	REAL,INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind) ! extinction coefficient in the visible
57	REAL,INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind)
58	REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth
59	! BENJAMIN MODIFS
60	real,intent(in) :: cloudfrac(ngrid,nlayer) ! cloud fraction
61	real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
62	logical,intent(in) :: clearsky
63
64	real aerosol0, obs_tau_col_aurora, pm, pente_cloud
65
66	INTEGER l,ig,iq,iaer
67
68	LOGICAL,SAVE :: firstcall=.true.
69	!$OMP THREADPRIVATE(firstcall)
70	REAL CBRT
71	EXTERNAL CBRT
72
73	INTEGER,SAVE :: i_co2ice=0 ! co2 ice
74	INTEGER,SAVE :: i_h2oice=0 ! water ice
75	!$OMP THREADPRIVATE(i_co2ice,i_h2oice)
76	CHARACTER(LEN=20) :: tracername ! to temporarily store text
77
78	! for fixed dust profiles
79	real topdust, expfactor, zp
80	REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling
81	REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling
82
83	real CLFtot
84
85	! identify tracers
86	IF (firstcall) THEN
87
88	write(,) "Tracers found in aeropacity:"
89	do iq=1,nq
90	tracername=noms(iq)
91	if (tracername.eq."co2_ice") then
92	i_co2ice=iq
93	write(,) "i_co2ice=",i_co2ice
94
95	endif
96	if (tracername.eq."h2o_ice") then
97	i_h2oice=iq
98	write(,) "i_h2oice=",i_h2oice
99	endif
100	enddo
101
102	if (noaero) then
103	print*, "No active aerosols found in aeropacity"
104	else
105	print*, "If you would like to use aerosols, make sure any old"
106	print*, "start files are updated in newstart using the option"
107	print*, "q=0"
108	write(,) "Active aerosols found in aeropacity:"
109	endif
110
111	if ((iaero_co2.ne.0).and.(.not.noaero)) then
112	print*, 'iaero_co2= ',iaero_co2
113	endif
114	if (iaero_h2o.ne.0) then
115	print*,'iaero_h2o= ',iaero_h2o
116	endif
117	if (iaero_dust.ne.0) then
118	print*,'iaero_dust= ',iaero_dust
119	endif
120	if (iaero_h2so4.ne.0) then
121	print*,'iaero_h2so4= ',iaero_h2so4
122	endif
123	if (iaero_back2lay.ne.0) then
124	print*,'iaero_back2lay= ',iaero_back2lay
125	endif
126	if (iaero_nh3.ne.0) then
127	print*,'iaero_nh3= ',iaero_nh3
128	endif
129	if (iaero_aurora.ne.0) then
130	print*,'iaero_aurora= ',iaero_aurora
131	endif
132
133	firstcall=.false.
134	ENDIF ! of IF (firstcall)
135
136
137	! ---------------------------------------------------------
138	!==================================================================
139	! CO2 ice aerosols
140	!==================================================================
141
142	if (iaero_co2.ne.0) then
143	iaer=iaero_co2
144	! 1. Initialization
145	aerosol(1:ngrid,1:nlayer,iaer)=0.0
146	! 2. Opacity calculation
147	if (noaero) then ! aerosol set to zero
148	aerosol(1:ngrid,1:nlayer,iaer)=0.0
149	elseif (aerofixco2.or.(i_co2ice.eq.0)) then ! CO2 ice cloud prescribed
150	aerosol(1:ngrid,1:nlayer,iaer)=1.e-9
151	!aerosol(1:ngrid,12,iaer)=4.0 ! single cloud layer option
152	else
153	DO ig=1, ngrid
154	DO l=1,nlayer-1 ! to stop the rad tran bug
155
156	aerosol0 = &
157	( 0.75 * QREFvis3d(ig,l,iaer) / &
158	( rho_co2 * reffrad(ig,l,iaer) ) ) * &
159	( pq(ig,l,i_co2ice) + 1.E-9 ) * &
160	( pplev(ig,l) - pplev(ig,l+1) ) / g
161	aerosol0 = max(aerosol0,1.e-9)
162	aerosol0 = min(aerosol0,L_TAUMAX)
163	aerosol(ig,l,iaer) = aerosol0
164	! aerosol(ig,l,iaer) = 0.0
165	! print*, aerosol(ig,l,iaer)
166	! using cloud fraction
167	! aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF))
168	! aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX)
169
170
171	ENDDO
172	ENDDO
173	end if ! if fixed or varying
174	end if ! if CO2 aerosols
175	!==================================================================
176	! Water ice / liquid
177	!==================================================================
178
179	if (iaero_h2o.ne.0) then
180	iaer=iaero_h2o
181	! 1. Initialization
182	aerosol(1:ngrid,1:nlayer,iaer)=0.0
183	! 2. Opacity calculation
184	if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
185	aerosol(1:ngrid,1:nlayer,iaer) =1.e-9
186
187	! put cloud at cloudlvl
188	if(kastprof.and.(cloudlvl.ne.0.0))then
189	ig=1
190	do l=1,nlayer
191	if(int(cloudlvl).eq.l)then
192	!if(cloudlvl.gt.(pplay(ig,l)/pplev(ig,1)))then
193	print*,'Inserting cloud at level ',l
194	!aerosol(ig,l,iaer)=10.0
195
196	rho_ice=920.0
197
198	! the Kasting approximation
199	aerosol(ig,l,iaer) = &
200	( 0.75 * QREFvis3d(ig,l,iaer) / &
201	( rho_ice * reffrad(ig,l,iaer) ) ) * &
202	!( pq(ig,l,i_h2oice) + 1.E-9 ) * &
203	( 4.0e-4 + 1.E-9 ) * &
204	( pplev(ig,l) - pplev(ig,l+1) ) / g
205
206
207	open(115,file='clouds.out',form='formatted')
208	write(115,*) l,aerosol(ig,l,iaer)
209	close(115)
210
211	return
212	endif
213	end do
214
215	call abort
216	endif
217
218	else
219
220	do ig=1, ngrid
221	do l=1,nlayer-1 ! to stop the rad tran bug
222
223	aerosol(ig,l,iaer) = & !modification by BC
224	( 0.75 * QREFvis3d(ig,l,iaer) / &
225	( rho_ice * reffrad(ig,l,iaer) ) ) * &
226	! pq(ig,l,i_h2oice) * & !JL I dropped the +1e-9 here to have the same
227	!( pplev(ig,l) - pplev(ig,l+1) ) / g ! opacity in the clearsky=true and the
228	! clear=false/pq=0 case
229	( pq(ig,l,i_h2oice) + 1.E-9 ) * & ! Doing this makes the code unstable, so I have restored it (RW)
230	( pplev(ig,l) - pplev(ig,l+1) ) / g
231
232	enddo
233	enddo
234
235	if(CLFvarying)then
236	call totalcloudfrac(ngrid,nlayer,nq,cloudfrac,totcloudfrac,pplev,pq,aerosol(1,1,iaer))
237	do ig=1, ngrid
238	do l=1,nlayer-1 ! to stop the rad tran bug
239	CLFtot = max(totcloudfrac(ig),0.01)
240	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
241	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
242	enddo
243	enddo
244	else
245	do ig=1, ngrid
246	do l=1,nlayer-1 ! to stop the rad tran bug
247	CLFtot = CLFfixval
248	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
249	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
250	enddo
251	enddo
252	end if!(CLFvarying)
253	endif !(aerofixed.or.(i_h2oice.eq.0).or.clearsky)
254
255	end if ! End if h2o aerosol
256
257	!==================================================================
258	! Dust
259	!==================================================================
260	if (iaero_dust.ne.0) then
261	iaer=iaero_dust
262	! 1. Initialization
263	aerosol(1:ngrid,1:nlayer,iaer)=0.0
264
265	topdust=30.0 ! km (used to be 10.0 km) LK
266
267	! 2. Opacity calculation
268
269	! expfactor=0.
270	DO l=1,nlayer-1
271	DO ig=1,ngrid
272	! Typical mixing ratio profile
273
274	zp=(pplev(ig,1)/pplay(ig,l))**(70./topdust)
275	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
276
277	! Vertical scaling function
278	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
279	*expfactor
280
281
282	ENDDO
283	ENDDO
284
285	! Rescaling each layer to reproduce the choosen (or assimilated)
286	! dust extinction opacity at visible reference wavelength, which
287	! is scaled to the surface pressure pplev(ig,1)
288
289	taudusttmp(1:ngrid)=0.
290	DO l=1,nlayer
291	DO ig=1,ngrid
292	taudusttmp(ig) = taudusttmp(ig) &
293	+ aerosol(ig,l,iaer)
294	ENDDO
295	ENDDO
296	DO l=1,nlayer-1
297	DO ig=1,ngrid
298	aerosol(ig,l,iaer) = max(1E-20, &
299	dusttau &
300	* pplev(ig,1) / pplev(ig,1) &
301	* aerosol(ig,l,iaer) &
302	/ taudusttmp(ig))
303
304	ENDDO
305	ENDDO
306	end if ! If dust aerosol
307
308	!==================================================================
309	! H2SO4
310	!==================================================================
311	! added by LK
312	if (iaero_h2so4.ne.0) then
313	iaer=iaero_h2so4
314
315	! 1. Initialization
316	aerosol(1:ngrid,1:nlayer,iaer)=0.0
317
318
319	! 2. Opacity calculation
320
321	! expfactor=0.
322	DO l=1,nlayer-1
323	DO ig=1,ngrid
324	! Typical mixing ratio profile
325
326	zp=(pplev(ig,1)/pplay(ig,l))**(70./30) !emulating topdust
327	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
328
329	! Vertical scaling function
330	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1))*expfactor
331
332	ENDDO
333	ENDDO
334	tauh2so4tmp(1:ngrid)=0.
335	DO l=1,nlayer
336	DO ig=1,ngrid
337	tauh2so4tmp(ig) = tauh2so4tmp(ig) + aerosol(ig,l,iaer)
338	ENDDO
339	ENDDO
340	DO l=1,nlayer-1
341	DO ig=1,ngrid
342	aerosol(ig,l,iaer) = max(1E-20, &
343	1 &
344	* pplev(ig,1) / pplev(ig,1) &
345	* aerosol(ig,l,iaer) &
346	/ tauh2so4tmp(ig))
347
348	ENDDO
349	ENDDO
350
351	! 1/700. is assuming a "sulfurtau" of 1
352	! Sulfur aerosol routine to be improved.
353	! aerosol0 = &
354	! ( 0.75 * QREFvis3d(ig,l,iaer) / &
355	! ( rho_h2so4 * reffrad(ig,l,iaer) ) ) * &
356	! ( pq(ig,l,i_h2so4) + 1.E-9 ) * &
357	! ( pplev(ig,l) - pplev(ig,l+1) ) / g
358	! aerosol0 = max(aerosol0,1.e-9)
359	! aerosol0 = min(aerosol0,L_TAUMAX)
360	! aerosol(ig,l,iaer) = aerosol0
361
362	! ENDDO
363	! ENDDO
364	end if
365
366
367	! ---------------------------------------------------------
368	!==================================================================
369	! Two-layer aerosols (unknown composition)
370	! S. Guerlet (2013)
371	!==================================================================
372
373	if (iaero_back2lay .ne.0) then
374	iaer=iaero_back2lay
375	! 1. Initialization
376	aerosol(1:ngrid,1:nlayer,iaer)=0.0
377	! 2. Opacity calculation
378	DO ig=1,ngrid
379	DO l=1,nlayer-1
380	aerosol(ig,l,iaer) = ( pplev(ig,l) - pplev(ig,l+1) )
381	!! 1. below tropospheric layer: no aerosols
382	IF (pplev(ig,l) .gt. pres_bottom_tropo) THEN
383	aerosol(ig,l,iaer) = 0.*aerosol(ig,l,iaer)
384	!! 2. tropo layer
385	ELSEIF (pplev(ig,l) .le. pres_bottom_tropo .and. pplev(ig,l) .ge. pres_top_tropo) THEN
386	aerosol(ig,l,iaer) = obs_tau_col_tropo*aerosol(ig,l,iaer)
387	!! 3. linear transition
388	ELSEIF (pplev(ig,l) .lt. pres_top_tropo .and. pplev(ig,l) .gt. pres_bottom_strato) THEN
389	expfactor=log(obs_tau_col_strato/obs_tau_col_tropo)/log(pres_bottom_strato/pres_top_tropo)
390	aerosol(ig,l,iaer)= obs_tau_col_tropo((pplev(ig,l)/pres_top_tropo)expfactor)aerosol(ig,l,iaer)/1.5
391	!! 4. strato layer
392	ELSEIF (pplev(ig,l) .le. pres_bottom_strato .and. pplev(ig,l) .gt. pres_top_strato) THEN
393	aerosol(ig,l,iaer)= obs_tau_col_strato*aerosol(ig,l,iaer)
394	!! 5. above strato layer: no aerosols
395	ELSEIF (pplev(ig,l) .lt. pres_top_strato) THEN
396	aerosol(ig,l,iaer) = 0.*aerosol(ig,l,iaer)
397	ENDIF
398	ENDDO
399	ENDDO
400
401	! 3. Re-normalize to observed total column
402	tau_col(:)=0.0
403	DO l=1,nlayer
404	DO ig=1,ngrid
405	tau_col(ig) = tau_col(ig) &
406	+ aerosol(ig,l,iaer)/(obs_tau_col_tropo+obs_tau_col_strato)
407	ENDDO
408	ENDDO
409
410	DO ig=1,ngrid
411	DO l=1,nlayer-1
412	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/tau_col(ig)
413	ENDDO
414	ENDDO
415
416
417	end if ! if Two-layer aerosols
418
419	!==================================================================
420	! Saturn/Jupiter ammonia cloud = thin cloud (scale height 0.2 hard coded...)
421	! S. Guerlet (2013)
422	!==================================================================
423
424	if (iaero_nh3 .ne.0) then
425	iaer=iaero_nh3
426	! 1. Initialization
427	aerosol(1:ngrid,1:nlayer,iaer)=0.D0
428	! 2. Opacity calculation
429	DO ig=1,ngrid
430
431	DO l=1,nlayer-1
432	!! 1. below cloud layer: no opacity
433
434	IF (pplev(ig,l) .gt. pres_nh3_cloud ) THEN
435	aerosol(ig,l,iaer) = 0.D0
436
437	ELSEIF (pplev(ig,l) .le. pres_nh3_cloud ) THEN
438	pente_cloud=5. !!(hard-coded, correspond to scale height 0.2)
439	aerosol(ig,l,iaer) = ((pplev(ig,l)/pres_nh3_cloud)*(pente_cloud))tau_nh3_cloud
440
441	ENDIF
442	ENDDO
443
444	END DO
445
446	! 3. Re-normalize to observed total column
447	tau_col(:)=0.0
448	DO l=1,nlayer
449	DO ig=1,ngrid
450	tau_col(ig) = tau_col(ig) &
451	+ aerosol(ig,l,iaer)/tau_nh3_cloud
452	ENDDO
453	ENDDO
454
455	DO ig=1,ngrid
456	DO l=1,nlayer-1
457	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/tau_col(ig)
458	ENDDO
459	ENDDO
460
461	end if ! if NH3 cloud
462	!==================================================================
463	! Jovian auroral aerosols (unknown composition) NON-GENERIC: vertical and meridional profile tuned to observations
464	! S. Guerlet (2015)
465	!==================================================================
466
467
468	if (iaero_aurora .ne.0) then
469	iaer=iaero_aurora
470	! 1. Initialization
471	aerosol(1:ngrid,1:nlayer,iaer)=0.D0
472	pm = 2000. !!case study: maxi aerosols at 20 hPa
473	! 2. Opacity calculation
474	DO ig=1,ngrid
475
476	!! Test Jupiter (based on Zhang et al 2013 observations, but a bit different), decembre 2015
477	DO l=1,nlayer
478	aerosol(ig,l,iaer) = (pplev(ig,l)/pm)*2 exp(-(pplev(ig,l)/pm)**2)
479	ENDDO
480	ENDDO
481
482	! 3. Meridional distribution, and re-normalize to observed total column
483	tau_col(:)=0.D0
484	DO ig=1,ngrid
485	!!Jupiter
486	!!Hem sud:
487	IF (latitude(ig)180.D0/pi .lt. -45.D0 .and. latitude(ig)180.D0/pi .gt. -70.) THEN
488	obs_tau_col_aurora= 10.D0*(-0.06D0latitude(ig)*180.D0/pi-3.4D0)
489	ELSEIF (latitude(ig)180.D0/pi .lt. -37.D0 .and. latitude(ig)180.D0/pi .ge. -45.) THEN
490	obs_tau_col_aurora= 10.D0*(-0.3D0latitude(ig)*180.D0/pi-14.3D0)
491	ELSEIF (latitude(ig)*180./pi .le. -70. ) THEN
492	obs_tau_col_aurora= 10*(0.0670.-3.4D0)
493	!!Hem Nord:
494	ELSEIF (latitude(ig)180.D0/pi .gt. 30.D0 .and. latitude(ig)180.D0/pi .lt. 70.) THEN
495	obs_tau_col_aurora= 10.D0*(0.03D0latitude(ig)*180.D0/pi-1.17D0)
496	ELSEIF (latitude(ig)180.D0/pi .gt. 22.D0 .and. latitude(ig)180.D0/pi .le. 30.) THEN
497	obs_tau_col_aurora= 10.D0*(0.3D0latitude(ig)*180.D0/pi-9.4D0)
498	ELSEIF (latitude(ig)*180.D0/pi .ge. 70.) THEN
499	obs_tau_col_aurora= 10*(0.0370.-1.17D0)
500	ELSEIF (latitude(ig)180.D0/pi .ge. -37. .and. latitude(ig)180.D0/pi .le. 22.) THEN
501	obs_tau_col_aurora = 0.001D0 !!Jupiter: mini pas a zero
502	ENDIF
503
504	DO l=1,nlayer
505	tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)/obs_tau_col_aurora
506	ENDDO
507	ENDDO
508
509	DO ig=1,ngrid
510	DO l=1,nlayer-1
511	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/tau_col(ig)
512	ENDDO
513	ENDDO
514
515
516	end if ! if Auroral aerosols
517
518
519	! --------------------------------------------------------------------------
520	! Column integrated visible optical depth in each point (used for diagnostic)
521
522	tau_col(:)=0.0
523	do iaer = 1, naerkind
524	do l=1,nlayer
525	do ig=1,ngrid
526	tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
527	end do
528	end do
529	end do
530
531	do ig=1,ngrid
532	do l=1,nlayer
533	do iaer = 1, naerkind
534	if(aerosol(ig,l,iaer).gt.1.e3)then
535	print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
536	print*,'at ig=',ig,', l=',l,', iaer=',iaer
537	print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
538	print*,'reffrad=',reffrad(ig,l,iaer)
539	endif
540	end do
541	end do
542	end do
543
544	do ig=1,ngrid
545	if(tau_col(ig).gt.1.e3)then
546	print*,'WARNING: tau_col=',tau_col(ig)
547	print*,'at ig=',ig
548	print*,'aerosol=',aerosol(ig,:,:)
549	print*,'QREFvis3d=',QREFvis3d(ig,:,:)
550	print*,'reffrad=',reffrad(ig,:,:)
551	endif
552	end do
553	return
554	end subroutine aeropacity
555

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

Download in other formats:

Original Format