Context Navigation

aeropacity.F90 @ 726

Last change on this file since 726 was 726, checked in by jleconte, 13 years ago

17/07/2012 == JL for LK

Generalization of aerosol scheme:
- any number of aerosols can be used and id numbers are determined consistently by the code. Aerosol order not important anymore.
- addition of a module with the id numbers for aerosols (aerosol_mod.F90).
- initialization of aerosols id numbers in iniaerosol.F90
- compile with -s x where x *must* be equal to the number of aerosols turned on in callphys.def (either by a flag or by dusttau>0 for dust). => may have to erase object files when compiling with s option for the first time.
For no aerosols, run with aeroco2=.true. and aerofixco2=.true (the default distribution for fixed co2

aerosols is 1.e-9; can be changed in aeropacity).

If starting from an old start file, recreate start file with the q=0 option in newstart.e.
update callphys.def with aeroXXX and aerofixXXX options (only XXX=co2,h2o supported for

now). Dust is activated by setting dusttau>0. See the early mars case in deftank.

To add other aerosols, see Laura Kerber.

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

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

Download in other formats:

Original Format