Context Navigation

aeropacity.F90 @ 255

Last change on this file since 255 was 253, checked in by emillour, 14 years ago

Generic GCM

Massive update to version 0.7

EM+RW

File size: 9.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 : 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	else
189	do ig=1, ngrid
190	do l=1,nlayer-1 ! to stop the rad tran bug
191
192	if(CLFvarying)then
193	CLF1 = max(cloudfrac(ig,l),0.01)
194	CLFtot = max(totcloudfrac(ig),0.01)
195	CLF2 = CLF1/CLFtot
196	CLF2 = min(1.,CLF2)
197	else
198	CLF1=1.0
199	CLF2=CLFfixval
200	endif
201
202	aerosol0 = &
203	( 0.75 * QREFvis3d(ig,l,iaer) / &
204	( rho_ice * reffrad(ig,l,iaer) ) ) * &
205	( pq(ig,l,i_h2oice) + 1.E-9 ) * &
206	( pplev(ig,l) - pplev(ig,l+1) ) / g / &
207	CLF1
208
209	aerosol(ig,l,iaer) = -log(1 - CLF2 + CLF2*exp(-aerosol0))
210	! why no division in exponential?
211	! because its already done in aerosol0
212
213	aerosol(ig,l,iaer) = max(aerosol(ig,l,iaer),1.e-9)
214	aerosol(ig,l,iaer) = min(aerosol(ig,l,iaer),aerosol0)
215
216	enddo
217	enddo
218	end if
219
220
221	!==================================================================
222	CASE(3) aerkind ! Dust (iaer=3)
223	!==================================================================
224
225	! 1. Initialization
226	aerosol(:,:,iaer)=0.0
227
228	topdust=10.0 ! km
229
230	print*,'WARNING, dust is experimental for Early Mars'
231
232	! 2. Opacity calculation
233
234	do l=1,nlayer
235	do ig=1,ngrid-1 ! to stop the rad tran bug
236	zp=700./pplay(ig,l)
237	aerosol(ig,l,1)=(dusttau/700.)*(pplev(ig,l)-pplev(ig,l+1)) &
238	max( exp(.03(1.-max(zp,1.))) , 1.E-3 ) &
239	*QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer)
240	! takes into account particle variation
241	enddo
242	enddo
243
244	!==================================================================
245	END SELECT aerkind
246	ENDDO ! iaer (loop on aerosol kind)
247
248
249	! --------------------------------------------------------------------------
250	! Column integrated visible optical depth in each point (used for diagnostic)
251
252	tau_col(:)=0.0
253	do iaer = 1, naerkind
254	do l=1,nlayer
255	do ig=1,ngrid
256	tau_col(ig) = tau_col(ig) + aerosol(ig,l,iaer)
257	end do
258	end do
259	end do
260
261	do ig=1, ngrid
262	do l=1,nlayer
263	do iaer = 1, naerkind
264	if(aerosol(ig,l,iaer).gt.1.e3)then
265	print*,'WARNING: aerosol=',aerosol(ig,l,iaer)
266	print*,'at ig=',ig,', l=',l,', iaer=',iaer
267	print*,'QREFvis3d=',QREFvis3d(ig,l,iaer)
268	print*,'reffrad=',reffrad(ig,l,iaer)
269	endif
270	end do
271	end do
272	end do
273
274	do ig=1, ngrid
275	if(tau_col(ig).gt.1.e3)then
276	print*,'WARNING: tau_col=',tau_col(ig)
277	print*,'at ig=',ig
278	print*,'aerosol=',aerosol(ig,:,:)
279	print*,'QREFvis3d=',QREFvis3d(ig,:,:)
280	print*,'reffrad=',reffrad(ig,:,:)
281	endif
282	end do
283
284	return
285	end subroutine aeropacity
286

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

Download in other formats:

Original Format