Context Navigation

aeropacity.F90 @ 1644

Last change on this file since 1644 was 1542, checked in by emillour, 9 years ago

Generic GCM:

fix for 1D in writediagfi to enable writing at "ecritphy" rate.
move iniprint.h to "misc"
Some code cleanup in anticipation of future updates:
- changed variable names in comgeomphy.F90: give them more explicit names: rlond => longitude , rlatd => latitude, airephy => cell_area, cuphy => dx , cvphy => dy
- removed long(), lati() and area() from comgeomfi_h.F90, use longitude(), latitude() and cell_are() from comgeomphy.F90 instead

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

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

Download in other formats:

Original Format