Context Navigation

aeropacity.F90 @ 601

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

13/02/2012 == JL + AS

All outputs are now in netCDF format. Even in 1D (No more G1D)
Clean up of the call to callcorrk when CLFvarying=true
Corrects a bug in writediagspecIR/VI. Output are now in W/m2/cm-1 as a function of the wavenumber in cm-1
Enable writediagspecIR/V to work in the CLFvarying=true case (output now done in Physiq after writediagfi)
Add a simple treatment for the supersaturation of CO2 (see forget et al 2012)
corrects a small bug when no clouds are present in aeropacity

File size: 10.5 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 : naerkind, L_TAUMAX
5
6	implicit none
7
8	!==================================================================
9	!
10	! Purpose
11	! -------
12	! Compute aerosol optical depth in each gridbox.
13	!
14	! Authors
15	! -------
16	! F. Forget
17	! F. Montmessin (water ice scheme)
18	! update J.-B. Madeleine (2008)
19	! dust removal, simplification by Robin Wordsworth (2009)
20	!
21	! Input
22	! -----
23	! ngrid Number of horizontal gridpoints
24	! nlayer Number of layers
25	! nq Number of tracers
26	! pplev Pressure (Pa) at each layer boundary
27	! pq Aerosol mixing ratio
28	! reffrad(ngrid,nlayer,naerkind) Aerosol effective radius
29	! QREFvis3d(ngridmx,nlayermx,naerkind) \ 3d extinction coefficients
30	! QREFir3d(ngridmx,nlayermx,naerkind) / at reference wavelengths
31	!
32	! Output
33	! ------
34	! aerosol Aerosol optical depth in layer l, grid point ig
35	! tau_col Total column optical depth at grid point ig
36	!
37	!=======================================================================
38
39	#include "dimensions.h"
40	#include "dimphys.h"
41	#include "callkeys.h"
42	#include "comcstfi.h"
43	#include "comgeomfi.h"
44	#include "tracer.h"
45
46	INTEGER ngrid,nlayer,nq
47
48	REAL pplay(ngrid,nlayer)
49	REAL pplev(ngrid,nlayer+1)
50	REAL pq(ngrid,nlayer,nq)
51	REAL aerosol(ngrid,nlayer,naerkind)
52	REAL reffrad(ngrid,nlayer,naerkind)
53	REAL QREFvis3d(ngridmx,nlayermx,naerkind)
54	REAL QREFir3d(ngridmx,nlayermx,naerkind)
55
56	REAL tau_col(ngrid)
57	! REAL tauref(ngrid), tau_col(ngrid)
58
59	real cloudfrac(ngridmx,nlayermx)
60	real aerosol0
61
62	INTEGER l,ig,iq,iaer
63
64	LOGICAL firstcall
65	DATA firstcall/.true./
66	SAVE firstcall
67
68	REAL CBRT
69	EXTERNAL CBRT
70
71	INTEGER,SAVE :: i_co2ice=0 ! co2 ice
72	INTEGER,SAVE :: i_h2oice=0 ! water ice
73	CHARACTER(LEN=20) :: tracername ! to temporarily store text
74
75	! for fixed dust profiles
76	real topdust, expfactor, zp
77	REAL taudusttmp(ngridmx) ! Temporary dust opacity used before scaling
78
79	! BENJAMIN MODIFS
80	real CLFtot,CLF1,CLF2
81	real totcloudfrac(ngridmx)
82	logical clearsky
83
84	! identify tracers
85	IF (firstcall) THEN
86
87	! are these tests of any real use ?
88	IF(ngrid.NE.ngridmx) THEN
89	PRINT*,'STOP in aeropacity!'
90	PRINT*,'problem of dimensions:'
91	PRINT*,'ngrid =',ngrid
92	PRINT*,'ngridmx =',ngridmx
93	STOP
94	ENDIF
95
96	if (nq.gt.nqmx) then
97	write(,) 'STOP in aeropacity: (nq .gt. nqmx)!'
98	write(,) 'nq=',nq,' nqmx=',nqmx
99	stop
100	endif
101
102	do iq=1,nqmx
103	tracername=noms(iq)
104	if (tracername.eq."co2_ice") then
105	i_co2ice=iq
106	endif
107	if (tracername.eq."h2o_ice") then
108	i_h2oice=iq
109	endif
110	enddo
111
112	write(,) "aeropacity: i_co2ice=",i_co2ice
113	write(,) " i_h2oice=",i_h2oice
114
115	if(watercond.and.(.not.aerofixed))then
116	if(naerkind.lt.2)then
117	print*,'Cannot have active H2O clouds with naerkind less than 2!'
118	call abort
119	endif
120	endif
121
122	if(dusttau.gt.0.0)then
123	if(naerkind.lt.3)then
124	print*,'Cannot have active dust with naerkind less than 3!'
125	call abort
126	endif
127	endif
128
129	firstcall=.false.
130	ENDIF ! of IF (firstcall)
131
132
133	DO iaer = 1, naerkind ! Loop on aerosol kind (NOT tracer)
134	! ---------------------------------------------------------
135	aerkind: SELECT CASE (iaer)
136	!==================================================================
137	CASE(1) aerkind ! CO2 ice (iaer=1)
138	!==================================================================
139
140	! 1. Initialization
141	aerosol(:,:,iaer)=0.0
142
143	! 2. Opacity calculation
144	if (aerofixed.or.(i_co2ice.eq.0)) then ! CO2 ice cloud prescribed
145	do ig=1, ngrid
146	do l=1,nlayer
147	aerosol(ig,l,iaer)=1.e-9
148	end do
149	!aerosol(ig,12,iaer)=4.0 ! single cloud layer option
150	end do
151	else
152	DO ig=1, ngrid
153	DO l=1,nlayer-1 ! to stop the rad tran bug
154
155	aerosol0 = &
156	( 0.75 * QREFvis3d(ig,l,iaer) / &
157	( rho_co2 * reffrad(ig,l,iaer) ) ) * &
158	( pq(ig,l,i_co2ice) + 1.E-9 ) * &
159	( pplev(ig,l) - pplev(ig,l+1) ) / g
160	aerosol0 = max(aerosol0,1.e-9)
161	aerosol0 = min(aerosol0,L_TAUMAX)
162	aerosol(ig,l,iaer) = aerosol0
163	! aerosol(ig,l,iaer) = 0.0
164
165	! using cloud fraction
166	! aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF))
167	! aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),L_TAUMAX)
168
169
170	ENDDO
171	ENDDO
172	end if
173
174	!==================================================================
175	CASE(2) aerkind ! Water ice / liquid (iaer=2)
176	!==================================================================
177
178	! 1. Initialization
179	aerosol(:,:,iaer)=0.0
180
181	! 2. Opacity calculation
182	if (aerofixed.or.(i_h2oice.eq.0).or.clearsky) then
183	do ig=1,ngrid ! to stop the rad tran bug
184	do l=1,nlayer
185	aerosol(ig,l,iaer) =1.e-9
186	end do
187	end do
188
189	! put cloud at cloudlvl
190	if(kastprof.and.(cloudlvl.ne.0.0))then
191	ig=1
192	do l=1,nlayer
193	if(int(cloudlvl).eq.l)then
194	!if(cloudlvl.gt.(pplay(ig,l)/pplev(ig,1)))then
195	print*,'Inserting cloud at level ',l
196	!aerosol(ig,l,iaer)=10.0
197
198	rho_ice=920.0
199
200	! the Kasting approximation
201	aerosol(ig,l,iaer) = &
202	( 0.75 * QREFvis3d(ig,l,iaer) / &
203	( rho_ice * reffrad(ig,l,iaer) ) ) * &
204	!( pq(ig,l,i_h2oice) + 1.E-9 ) * &
205	( 4.0e-4 + 1.E-9 ) * &
206	( pplev(ig,l) - pplev(ig,l+1) ) / g
207
208
209	open(115,file='clouds.out',form='formatted')
210	write(115,*) l,aerosol(ig,l,iaer)
211	close(115)
212
213	return
214	endif
215	end do
216
217	call abort
218	endif
219
220	else
221	do ig=1, ngrid
222	do l=1,nlayer-1 ! to stop the rad tran bug
223
224	if(CLFvarying)then
225	CLF1 = max(cloudfrac(ig,l),0.01)
226	CLFtot = max(totcloudfrac(ig),0.01)
227	CLF2 = CLF1/CLFtot
228	CLF2 = min(1.,CLF2)
229	else
230	CLF1=1.0
231	CLF2=CLFfixval
232	endif
233
234	aerosol0 = &
235	( 0.75 * QREFvis3d(ig,l,iaer) / &
236	( rho_ice * reffrad(ig,l,iaer) ) ) * &
237	pq(ig,l,i_h2oice) * & !JL I dropped the +1e-9 here to have the same
238	( pplev(ig,l) - pplev(ig,l+1) ) / g / & ! opacity in the clearsky=true and the
239	CLF1 ! clear=false/pq=0 case
240
241	aerosol(ig,l,iaer) = -log(1 - CLF2 + CLF2*exp(-aerosol0))
242	! why no division in exponential?
243	! because its already done in aerosol0
244
245	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
246	aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),aerosol0)
247
248	enddo
249	enddo
250	end if
251
252
253	!==================================================================
254	CASE(3) aerkind ! Dust (iaer=3)
255	!==================================================================
256
257	! 1. Initialization
258	aerosol(:,:,iaer)=0.0
259
260	topdust=10.0 ! km
261
262	print*,'WARNING, dust is experimental for Early Mars'
263
264	! 2. Opacity calculation
265
266	do ig=1,ngrid
267	do l=1,nlayer
268	zp=700./pplay(ig,l)
269	aerosol(ig,l,iaer)=(dusttau/700.)*(pplev(ig,l)-pplev(ig,l+1)) &
270	max( exp(.03(1.-max(zp,1.))) , 1.E-3 ) &
271	*QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer)
272	! takes into account particle variation
273	enddo
274	enddo
275
276	!==================================================================
277	END SELECT aerkind
278	ENDDO ! iaer (loop on aerosol kind)
279
280
281	! --------------------------------------------------------------------------
282	! Column integrated visible optical depth in each point (used for diagnostic)
283
284	tau_col(:)=0.0
285	do iaer = 1, naerkind
286	do l=1,nlayer
287	do ig=1,ngrid
288	tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
289	end do
290	end do
291	end do
292
293	do ig=1, ngrid
294	do l=1,nlayer
295	do iaer = 1, naerkind
296	if(aerosol(ig,l,iaer).gt.1.e3)then
297	print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
298	print*,'at ig=',ig,', l=',l,', iaer=',iaer
299	print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
300	print*,'reffrad=',reffrad(ig,l,iaer)
301	endif
302	end do
303	end do
304	end do
305
306	do ig=1, ngrid
307	if(tau_col(ig).gt.1.e3)then
308	print*,'WARNING: tau_col=',tau_col(ig)
309	print*,'at ig=',ig
310	print*,'aerosol=',aerosol(ig,:,:)
311	print*,'QREFvis3d=',QREFvis3d(ig,:,:)
312	print*,'reffrad=',reffrad(ig,:,:)
313	endif
314	end do
315
316	return
317	end subroutine aeropacity
318

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

Download in other formats:

Original Format