Context Navigation

aeropacity.F90 @ 724

Last change on this file since 724 was 716, checked in by rwordsworth, 13 years ago

Mainly updates to radiative transfer and gas management scheme.
Most CIA data now read from standard HITRAN datafiles. For the H2O
continuum, two options have been added: the standard CKD continuum,
and the empirical formula in PPC (Pierrehumbert 2010). Use the toggle
'H2Ocont_simple' in callphys.def to choose.

Note to Martians: I've changed the default values of 'sedimentation' and
'co2cond' in inifis to false. Both these are defined in the standard deftank
callphys.def file, so there shouldn't be any compatibility problems.

File size: 10.7 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	( pq(ig,l,i_h2oice) + 1.E-9 ) * & ! Doing this makes the code unstable, so I have restored it (RW)
241	( pplev(ig,l) - pplev(ig,l+1) ) / g / &
242	CLF1
243
244	aerosol(ig,l,iaer) = -log(1 - CLF2 + CLF2*exp(-aerosol0))
245	! why no division in exponential?
246	! because its already done in aerosol0
247
248	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
249	aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),aerosol0)
250
251	enddo
252	enddo
253	end if
254
255
256	!==================================================================
257	CASE(3) aerkind ! Dust (iaer=3)
258	!==================================================================
259
260	! 1. Initialization
261	aerosol(:,:,iaer)=0.0
262
263	topdust=10.0 ! km
264
265	print*,'WARNING, dust is experimental for Early Mars'
266
267	! 2. Opacity calculation
268
269	do ig=1,ngrid
270	do l=1,nlayer
271	zp=700./pplay(ig,l)
272	aerosol(ig,l,iaer)=(dusttau/700.)*(pplev(ig,l)-pplev(ig,l+1)) &
273	max( exp(.03(1.-max(zp,1.))) , 1.E-3 ) &
274	*QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer)
275	! takes into account particle variation
276	enddo
277	enddo
278
279	!==================================================================
280	END SELECT aerkind
281	ENDDO ! iaer (loop on aerosol kind)
282
283
284	! --------------------------------------------------------------------------
285	! Column integrated visible optical depth in each point (used for diagnostic)
286
287	tau_col(:)=0.0
288	do iaer = 1, naerkind
289	do l=1,nlayer
290	do ig=1,ngrid
291	tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
292	end do
293	end do
294	end do
295
296	do ig=1, ngrid
297	do l=1,nlayer
298	do iaer = 1, naerkind
299	if(aerosol(ig,l,iaer).gt.1.e3)then
300	print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
301	print*,'at ig=',ig,', l=',l,', iaer=',iaer
302	print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
303	print*,'reffrad=',reffrad(ig,l,iaer)
304	endif
305	end do
306	end do
307	end do
308
309	do ig=1, ngrid
310	if(tau_col(ig).gt.1.e3)then
311	print*,'WARNING: tau_col=',tau_col(ig)
312	print*,'at ig=',ig
313	print*,'aerosol=',aerosol(ig,:,:)
314	print*,'QREFvis3d=',QREFvis3d(ig,:,:)
315	print*,'reffrad=',reffrad(ig,:,:)
316	endif
317	end do
318
319	return
320	end subroutine aeropacity
321

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

Download in other formats:

Original Format