Context Navigation

aeropacity.F90 @ 306

Last change on this file since 306 was 305, checked in by rwordsworth, 14 years ago
Several new files added as part of the climate evolution model (main program kcm.F90). Some general cleanup in physiq.F90 and callcorrk.F90. Bugs in dust radiative transfer and H2 Rayleigh scattering corrected.
File size: 10.3 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) + 1.E-9 ) * &
238	( pplev(ig,l) - pplev(ig,l+1) ) / g / &
239	CLF1
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,1)=(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