Context Navigation

aeropacity.F90 @ 764

Last change on this file since 764 was 741, checked in by jleconte, 13 years ago
Model now works without aerosols at all. You still need to compile it with -s 1 option. LK12
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
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 (noaero) then
119	print*, "No active aerosols found in aeropacity"
120	else
121	print*, "If you would like to use aerosols, make sure any old"
122	print*, "start files are updated in newstart using the option"
123	print*, "q=0"
124	write(,) "Active aerosols found in aeropacity:"
125	endif
126
127	if ((iaero_co2.ne.0).and.(.not.noaero)) then
128	print*, 'iaero_co2= ',iaero_co2
129	endif
130	if (iaero_h2o.ne.0) then
131	print*,'iaero_h2o= ',iaero_h2o
132	endif
133	if (iaero_dust.ne.0) then
134	print*,'iaero_dust= ',iaero_dust
135	endif
136	if (iaero_h2so4.ne.0) then
137	print*,'iaero_h2so4= ',iaero_h2so4
138	endif
139
140	firstcall=.false.
141	ENDIF ! of IF (firstcall)
142
143
144	! ---------------------------------------------------------
145	!==================================================================
146	! CO2 ice aerosols
147	!==================================================================
148
149	if (iaero_co2.ne.0) then
150	iaer=iaero_co2
151	! 1. Initialization
152	aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
153	! 2. Opacity calculation
154	if (noaero) then ! aerosol set to a very low value
155	aerosol(1:ngridmx,1:nlayermx,iaer)=1.e-9
156	elseif (aerofixco2.or.(i_co2ice.eq.0)) then ! CO2 ice cloud prescribed
157	aerosol(1:ngridmx,1:nlayermx,iaer)=1.e-9
158	!aerosol(1:ngridmx,12,iaer)=4.0 ! single cloud layer option
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(1:ngridmx,1:nlayermx,iaer)=0.0
190	! 2. Opacity calculation
191	if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
192	aerosol(1:ngridmx,1:nlayermx,iaer) =1.e-9
193
194	! put cloud at cloudlvl
195	if(kastprof.and.(cloudlvl.ne.0.0))then
196	ig=1
197	do l=1,nlayer
198	if(int(cloudlvl).eq.l)then
199	!if(cloudlvl.gt.(pplay(ig,l)/pplev(ig,1)))then
200	print*,'Inserting cloud at level ',l
201	!aerosol(ig,l,iaer)=10.0
202
203	rho_ice=920.0
204
205	! the Kasting approximation
206	aerosol(ig,l,iaer) = &
207	( 0.75 * QREFvis3d(ig,l,iaer) / &
208	( rho_ice * reffrad(ig,l,iaer) ) ) * &
209	!( pq(ig,l,i_h2oice) + 1.E-9 ) * &
210	( 4.0e-4 + 1.E-9 ) * &
211	( pplev(ig,l) - pplev(ig,l+1) ) / g
212
213
214	open(115,file='clouds.out',form='formatted')
215	write(115,*) l,aerosol(ig,l,iaer)
216	close(115)
217
218	return
219	endif
220	end do
221
222	call abort
223	endif
224
225	else
226
227	do ig=1, ngrid
228	do l=1,nlayer-1 ! to stop the rad tran bug
229
230	aerosol(ig,l,iaer) = & !modification by BC
231	( 0.75 * QREFvis3d(ig,l,iaer) / &
232	( rho_ice * reffrad(ig,l,iaer) ) ) * &
233	! pq(ig,l,i_h2oice) * & !JL I dropped the +1e-9 here to have the same
234	!( pplev(ig,l) - pplev(ig,l+1) ) / g ! opacity in the clearsky=true and the
235	! clear=false/pq=0 case
236	( pq(ig,l,i_h2oice) + 1.E-9 ) * & ! Doing this makes the code unstable, so I have restored it (RW)
237	( pplev(ig,l) - pplev(ig,l+1) ) / g
238
239	enddo
240	enddo
241
242	if(CLFvarying)then
243	call totalcloudfrac(cloudfrac,totcloudfrac,pplev,pq,aerosol(1:ngridmx,1:nlayermx,iaer))
244	do ig=1, ngrid
245	do l=1,nlayer-1 ! to stop the rad tran bug
246	CLFtot = max(totcloudfrac(ig),0.01)
247	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
248	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
249	enddo
250	enddo
251	else
252	do ig=1, ngrid
253	do l=1,nlayer-1 ! to stop the rad tran bug
254	CLFtot = CLFfixval
255	aerosol(ig,l,iaer)=aerosol(ig,l,iaer)/CLFtot
256	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
257	enddo
258	enddo
259	end if!(CLFvarying)
260	endif !(aerofixed.or.(i_h2oice.eq.0).or.clearsky)
261
262	end if ! End if h2o aerosol
263
264	!==================================================================
265	! Dust
266	!==================================================================
267	if (iaero_dust.ne.0) then
268	iaer=iaero_dust
269	! 1. Initialization
270	aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
271
272	topdust=30.0 ! km (used to be 10.0 km) LK
273
274	! 2. Opacity calculation
275
276	! expfactor=0.
277	DO l=1,nlayer-1
278	DO ig=1,ngrid
279	! Typical mixing ratio profile
280
281	zp=(preff/pplay(ig,l))**(70./topdust)
282	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
283
284	! Vertical scaling function
285	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
286	*expfactor
287
288
289	ENDDO
290	ENDDO
291
292	! Rescaling each layer to reproduce the choosen (or assimilated)
293	! dust extinction opacity at visible reference wavelength, which
294	! is scaled to the "preff" reference surface pressure available in comvert.h
295	! and stored in startfi.nc
296
297	taudusttmp(1:ngrid)=0.
298	DO l=1,nlayer
299	DO ig=1,ngrid
300	taudusttmp(ig) = taudusttmp(ig) &
301	+ aerosol(ig,l,iaer)
302	ENDDO
303	ENDDO
304	DO l=1,nlayer-1
305	DO ig=1,ngrid
306	aerosol(ig,l,iaer) = max(1E-20, &
307	dusttau &
308	* pplev(ig,1) / preff &
309	* aerosol(ig,l,iaer) &
310	/ taudusttmp(ig))
311
312	ENDDO
313	ENDDO
314	end if ! If dust aerosol
315
316	!==================================================================
317	! H2SO4
318	!==================================================================
319	! added by LK
320	if (iaero_h2so4.ne.0) then
321	iaer=iaero_h2so4
322
323	! 1. Initialization
324	aerosol(1:ngridmx,1:nlayermx,iaer)=0.0
325
326
327	! 2. Opacity calculation
328
329	! expfactor=0.
330	DO l=1,nlayer-1
331	DO ig=1,ngrid
332	! Typical mixing ratio profile
333
334	zp=(preff/pplay(ig,l))**(70./30) !emulating topdust
335	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
336
337	! Vertical scaling function
338	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1))*expfactor
339
340	ENDDO
341	ENDDO
342	tauh2so4tmp(1:ngrid)=0.
343	DO l=1,nlayer
344	DO ig=1,ngrid
345	tauh2so4tmp(ig) = tauh2so4tmp(ig) + aerosol(ig,l,iaer)
346	ENDDO
347	ENDDO
348	DO l=1,nlayer-1
349	DO ig=1,ngrid
350	aerosol(ig,l,iaer) = max(1E-20, &
351	1 &
352	* pplev(ig,1) / preff &
353	* aerosol(ig,l,iaer) &
354	/ tauh2so4tmp(ig))
355
356	ENDDO
357	ENDDO
358
359	! 1/700. is assuming a "sulfurtau" of 1
360	! Sulfur aerosol routine to be improved.
361	! aerosol0 = &
362	! ( 0.75 * QREFvis3d(ig,l,iaer) / &
363	! ( rho_h2so4 * reffrad(ig,l,iaer) ) ) * &
364	! ( pq(ig,l,i_h2so4) + 1.E-9 ) * &
365	! ( pplev(ig,l) - pplev(ig,l+1) ) / g
366	! aerosol0 = max(aerosol0,1.e-9)
367	! aerosol0 = min(aerosol0,L_TAUMAX)
368	! aerosol(ig,l,iaer) = aerosol0
369
370	! ENDDO
371	! ENDDO
372	end if
373
374
375
376	! --------------------------------------------------------------------------
377	! Column integrated visible optical depth in each point (used for diagnostic)
378
379	tau_col(:)=0.0
380	do iaer = 1, naerkind
381	do l=1,nlayer
382	do ig=1,ngrid
383	tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
384	end do
385	end do
386	end do
387
388	do ig=1, ngrid
389	do l=1,nlayer
390	do iaer = 1, naerkind
391	if(aerosol(ig,l,iaer).gt.1.e3)then
392	print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
393	print*,'at ig=',ig,', l=',l,', iaer=',iaer
394	print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
395	print*,'reffrad=',reffrad(ig,l,iaer)
396	endif
397	end do
398	end do
399	end do
400
401	do ig=1, ngrid
402	if(tau_col(ig).gt.1.e3)then
403	print*,'WARNING: tau_col=',tau_col(ig)
404	print*,'at ig=',ig
405	print*,'aerosol=',aerosol(ig,:,:)
406	print*,'QREFvis3d=',QREFvis3d(ig,:,:)
407	print*,'reffrad=',reffrad(ig,:,:)
408	endif
409	end do
410	return
411	end subroutine aeropacity
412

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

Download in other formats:

Original Format