Context Navigation

aeropacity.F90 @ 867

Last change on this file since 867 was 858, checked in by emillour, 12 years ago

Generic GCM:

Fixed sedimentation issue: ensure in callsedim that the correct radii are provided to newsedim and also that the updated temperature and tracer mixing ratios are used to compute sedimentation.
Updated the way aerosol radii are considered and used; routines in radii_mod (h2o_reffrad, co2_reffrad, etc.) only handle a single aerosol. The idea here is that these can be called from anywhere and that the caller doesn't need to have the full (naerkind size) array of aerosol radii.
cleanup (addition of intent(..) to routine arguments) in various routines

File size: 13.9 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 comgeomfi_h
7	USE tracer_h, only: noms,rho_co2,rho_ice
8
9	implicit none
10
11	!==================================================================
12	!
13	! Purpose
14	! -------
15	! Compute aerosol optical depth in each gridbox.
16	!
17	! Authors
18	! -------
19	! F. Forget
20	! F. Montmessin (water ice scheme)
21	! update J.-B. Madeleine (2008)
22	! dust removal, simplification by Robin Wordsworth (2009)
23	!
24	! Input
25	! -----
26	! ngrid Number of horizontal gridpoints
27	! nlayer Number of layers
28	! nq Number of tracers
29	! pplev Pressure (Pa) at each layer boundary
30	! pq Aerosol mixing ratio
31	! reffrad(ngrid,nlayer,naerkind) Aerosol effective radius
32	! QREFvis3d(ngrid,nlayermx,naerkind) \ 3d extinction coefficients
33	! QREFir3d(ngrid,nlayermx,naerkind) / at reference wavelengths
34	!
35	! Output
36	! ------
37	! aerosol Aerosol optical depth in layer l, grid point ig
38	! tau_col Total column optical depth at grid point ig
39	!
40	!=======================================================================
41
42	#include "dimensions.h"
43	#include "dimphys.h"
44	#include "callkeys.h"
45	#include "comcstfi.h"
46	#include "comvert.h"
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,nlayermx,naerkind) ! extinction coefficient in the visible
57	REAL,INTENT(IN) :: QREFir3d(ngrid,nlayermx,naerkind)
58	REAL,INTENT(OUT):: tau_col(ngrid) !column integrated visible optical depth
59	! BENJAMIN MODIFS
60	real,intent(in) :: cloudfrac(ngrid,nlayermx) ! cloud fraction
61	real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
62	logical,intent(in) :: clearsky
63
64	real aerosol0
65
66	INTEGER l,ig,iq,iaer
67
68	LOGICAL,SAVE :: firstcall=.true.
69
70	REAL CBRT
71	EXTERNAL CBRT
72
73	INTEGER,SAVE :: i_co2ice=0 ! co2 ice
74	INTEGER,SAVE :: i_h2oice=0 ! water ice
75	CHARACTER(LEN=20) :: tracername ! to temporarily store text
76
77	! for fixed dust profiles
78	real topdust, expfactor, zp
79	REAL taudusttmp(ngrid) ! Temporary dust opacity used before scaling
80	REAL tauh2so4tmp(ngrid) ! Temporary h2so4 opacity used before scaling
81
82	real CLFtot
83
84	! identify tracers
85	IF (firstcall) THEN
86
87	write(,) "Tracers found in aeropacity:"
88	do iq=1,nq
89	tracername=noms(iq)
90	if (tracername.eq."co2_ice") then
91	i_co2ice=iq
92	write(,) "i_co2ice=",i_co2ice
93
94	endif
95	if (tracername.eq."h2o_ice") then
96	i_h2oice=iq
97	write(,) "i_h2oice=",i_h2oice
98	endif
99	enddo
100
101	if (noaero) then
102	print*, "No active aerosols found in aeropacity"
103	else
104	print*, "If you would like to use aerosols, make sure any old"
105	print*, "start files are updated in newstart using the option"
106	print*, "q=0"
107	write(,) "Active aerosols found in aeropacity:"
108	endif
109
110	if ((iaero_co2.ne.0).and.(.not.noaero)) then
111	print*, 'iaero_co2= ',iaero_co2
112	endif
113	if (iaero_h2o.ne.0) then
114	print*,'iaero_h2o= ',iaero_h2o
115	endif
116	if (iaero_dust.ne.0) then
117	print*,'iaero_dust= ',iaero_dust
118	endif
119	if (iaero_h2so4.ne.0) then
120	print*,'iaero_h2so4= ',iaero_h2so4
121	endif
122
123	firstcall=.false.
124	ENDIF ! of IF (firstcall)
125
126
127	! ---------------------------------------------------------
128	!==================================================================
129	! CO2 ice aerosols
130	!==================================================================
131
132	if (iaero_co2.ne.0) then
133	iaer=iaero_co2
134	! 1. Initialization
135	aerosol(1:ngrid,1:nlayermx,iaer)=0.0
136	! 2. Opacity calculation
137	if (noaero) then ! aerosol set to zero
138	aerosol(1:ngrid,1:nlayermx,iaer)=0.0
139	elseif (aerofixco2.or.(i_co2ice.eq.0)) then ! CO2 ice cloud prescribed
140	aerosol(1:ngrid,1:nlayermx,iaer)=1.e-9
141	!aerosol(1:ngrid,12,iaer)=4.0 ! single cloud layer option
142	else
143	DO ig=1, ngrid
144	DO l=1,nlayer-1 ! to stop the rad tran bug
145
146	aerosol0 = &
147	( 0.75 * QREFvis3d(ig,l,iaer) / &
148	( rho_co2 * reffrad(ig,l,iaer) ) ) * &
149	( pq(ig,l,i_co2ice) + 1.E-9 ) * &
150	( pplev(ig,l) - pplev(ig,l+1) ) / g
151	aerosol0 = max(aerosol0,1.e-9)
152	aerosol0 = min(aerosol0,L_TAUMAX)
153	aerosol(ig,l,iaer) = aerosol0
154	! aerosol(ig,l,iaer) = 0.0
155	! print*, aerosol(ig,l,iaer)
156	! using cloud fraction
157	! aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF))
158	! aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX)
159
160
161	ENDDO
162	ENDDO
163	end if ! if fixed or varying
164	end if ! if CO2 aerosols
165	!==================================================================
166	! Water ice / liquid
167	!==================================================================
168
169	if (iaero_h2o.ne.0) then
170	iaer=iaero_h2o
171	! 1. Initialization
172	aerosol(1:ngrid,1:nlayermx,iaer)=0.0
173	! 2. Opacity calculation
174	if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
175	aerosol(1:ngrid,1:nlayermx,iaer) =1.e-9
176
177	! put cloud at cloudlvl
178	if(kastprof.and.(cloudlvl.ne.0.0))then
179	ig=1
180	do l=1,nlayer
181	if(int(cloudlvl).eq.l)then
182	!if(cloudlvl.gt.(pplay(ig,l)/pplev(ig,1)))then
183	print*,'Inserting cloud at level ',l
184	!aerosol(ig,l,iaer)=10.0
185
186	rho_ice=920.0
187
188	! the Kasting approximation
189	aerosol(ig,l,iaer) = &
190	( 0.75 * QREFvis3d(ig,l,iaer) / &
191	( rho_ice * reffrad(ig,l,iaer) ) ) * &
192	!( pq(ig,l,i_h2oice) + 1.E-9 ) * &
193	( 4.0e-4 + 1.E-9 ) * &
194	( pplev(ig,l) - pplev(ig,l+1) ) / g
195
196
197	open(115,file='clouds.out',form='formatted')
198	write(115,*) l,aerosol(ig,l,iaer)
199	close(115)
200
201	return
202	endif
203	end do
204
205	call abort
206	endif
207
208	else
209
210	do ig=1, ngrid
211	do l=1,nlayer-1 ! to stop the rad tran bug
212
213	aerosol(ig,l,iaer) = & !modification by BC
214	( 0.75 * QREFvis3d(ig,l,iaer) / &
215	( rho_ice * reffrad(ig,l,iaer) ) ) * &
216	! pq(ig,l,i_h2oice) * & !JL I dropped the +1e-9 here to have the same
217	!( pplev(ig,l) - pplev(ig,l+1) ) / g ! opacity in the clearsky=true and the
218	! clear=false/pq=0 case
219	( pq(ig,l,i_h2oice) + 1.E-9 ) * & ! Doing this makes the code unstable, so I have restored it (RW)
220	( pplev(ig,l) - pplev(ig,l+1) ) / g
221
222	enddo
223	enddo
224
225	if(CLFvarying)then
226	call totalcloudfrac(ngrid,nq,cloudfrac,totcloudfrac,pplev,pq,aerosol(1:ngrid,1:nlayermx,iaer))
227	do ig=1, ngrid
228	do l=1,nlayer-1 ! to stop the rad tran bug
229	CLFtot = max(totcloudfrac(ig),0.01)
230	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
231	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
232	enddo
233	enddo
234	else
235	do ig=1, ngrid
236	do l=1,nlayer-1 ! to stop the rad tran bug
237	CLFtot = CLFfixval
238	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
239	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
240	enddo
241	enddo
242	end if!(CLFvarying)
243	endif !(aerofixed.or.(i_h2oice.eq.0).or.clearsky)
244
245	end if ! End if h2o aerosol
246
247	!==================================================================
248	! Dust
249	!==================================================================
250	if (iaero_dust.ne.0) then
251	iaer=iaero_dust
252	! 1. Initialization
253	aerosol(1:ngrid,1:nlayermx,iaer)=0.0
254
255	topdust=30.0 ! km (used to be 10.0 km) LK
256
257	! 2. Opacity calculation
258
259	! expfactor=0.
260	DO l=1,nlayer-1
261	DO ig=1,ngrid
262	! Typical mixing ratio profile
263
264	zp=(preff/pplay(ig,l))**(70./topdust)
265	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
266
267	! Vertical scaling function
268	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
269	*expfactor
270
271
272	ENDDO
273	ENDDO
274
275	! Rescaling each layer to reproduce the choosen (or assimilated)
276	! dust extinction opacity at visible reference wavelength, which
277	! is scaled to the "preff" reference surface pressure available in comvert.h
278	! and stored in startfi.nc
279
280	taudusttmp(1:ngrid)=0.
281	DO l=1,nlayer
282	DO ig=1,ngrid
283	taudusttmp(ig) = taudusttmp(ig) &
284	+ aerosol(ig,l,iaer)
285	ENDDO
286	ENDDO
287	DO l=1,nlayer-1
288	DO ig=1,ngrid
289	aerosol(ig,l,iaer) = max(1E-20, &
290	dusttau &
291	* pplev(ig,1) / preff &
292	* aerosol(ig,l,iaer) &
293	/ taudusttmp(ig))
294
295	ENDDO
296	ENDDO
297	end if ! If dust aerosol
298
299	!==================================================================
300	! H2SO4
301	!==================================================================
302	! added by LK
303	if (iaero_h2so4.ne.0) then
304	iaer=iaero_h2so4
305
306	! 1. Initialization
307	aerosol(1:ngrid,1:nlayermx,iaer)=0.0
308
309
310	! 2. Opacity calculation
311
312	! expfactor=0.
313	DO l=1,nlayer-1
314	DO ig=1,ngrid
315	! Typical mixing ratio profile
316
317	zp=(preff/pplay(ig,l))**(70./30) !emulating topdust
318	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
319
320	! Vertical scaling function
321	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1))*expfactor
322
323	ENDDO
324	ENDDO
325	tauh2so4tmp(1:ngrid)=0.
326	DO l=1,nlayer
327	DO ig=1,ngrid
328	tauh2so4tmp(ig) = tauh2so4tmp(ig) + aerosol(ig,l,iaer)
329	ENDDO
330	ENDDO
331	DO l=1,nlayer-1
332	DO ig=1,ngrid
333	aerosol(ig,l,iaer) = max(1E-20, &
334	1 &
335	* pplev(ig,1) / preff &
336	* aerosol(ig,l,iaer) &
337	/ tauh2so4tmp(ig))
338
339	ENDDO
340	ENDDO
341
342	! 1/700. is assuming a "sulfurtau" of 1
343	! Sulfur aerosol routine to be improved.
344	! aerosol0 = &
345	! ( 0.75 * QREFvis3d(ig,l,iaer) / &
346	! ( rho_h2so4 * reffrad(ig,l,iaer) ) ) * &
347	! ( pq(ig,l,i_h2so4) + 1.E-9 ) * &
348	! ( pplev(ig,l) - pplev(ig,l+1) ) / g
349	! aerosol0 = max(aerosol0,1.e-9)
350	! aerosol0 = min(aerosol0,L_TAUMAX)
351	! aerosol(ig,l,iaer) = aerosol0
352
353	! ENDDO
354	! ENDDO
355	end if
356
357
358
359	! --------------------------------------------------------------------------
360	! Column integrated visible optical depth in each point (used for diagnostic)
361
362	tau_col(:)=0.0
363	do iaer = 1, naerkind
364	do l=1,nlayer
365	do ig=1,ngrid
366	tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
367	end do
368	end do
369	end do
370
371	do ig=1,ngrid
372	do l=1,nlayer
373	do iaer = 1, naerkind
374	if(aerosol(ig,l,iaer).gt.1.e3)then
375	print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
376	print*,'at ig=',ig,', l=',l,', iaer=',iaer
377	print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
378	print*,'reffrad=',reffrad(ig,l,iaer)
379	endif
380	end do
381	end do
382	end do
383
384	do ig=1,ngrid
385	if(tau_col(ig).gt.1.e3)then
386	print*,'WARNING: tau_col=',tau_col(ig)
387	print*,'at ig=',ig
388	print*,'aerosol=',aerosol(ig,:,:)
389	print*,'QREFvis3d=',QREFvis3d(ig,:,:)
390	print*,'reffrad=',reffrad(ig,:,:)
391	endif
392	end do
393	return
394	end subroutine aeropacity
395

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

Download in other formats:

Original Format