Context Navigation

aeropacity.F90 @ 728

Last change on this file since 728 was 728, checked in by jleconte, 12 years ago

18/07/2012 == JL

New water cycle scheme:
- largescale now in F90. Robustness increased by i) including evap inside largescale ii) computing the

condensed water amount iteratively

same improvements in moistadj.
Water thermodynamical data and saturation curves centralized in module watercommn_h
- The saturation curves used are now Tetens formula as they are analyticaly inversible (Ts(P)-> Ps(T)).

New saturation curve yields very good agreement with the former one.

Saturation curves are now generalized for arbitrary water amount (not just q<<1)
The old watersat should be removed soon.

The effect of water vapor on total (surface) pressure can be taken into account by setting

mass_redistrib=.true. in callphys.def (routine mass_redistribution inspired from co2_condense in martian
model but with a different scheme as many routines evaporate/condense water vapor).

New cloud and precipitation scheme (JL + BC):
- The default recovery assumption for computing the total cloud fraction has been changed (total random gave too

large cloud fractions). See totalcloudfrac.F90 for details and to change this.

Totalcloudfraction now set the total cloud fraction to the fraction of the

optically thickest cloud and totalcloudfrac is thus called in aeropacity.

Only the total cloud fraction is used to compute optical depth in aeropacity (no more effective

optical depth with exponential formula).

4 precipitation schemes are now available (see rain.F90 for details). The choice can be made using precip_scheme

in callphys.def. Usage of the more physically based model of Boucher et al 95 (precip_scheme=4) is recommended.
default behavior is set to the former "simple scheme" (precip_scheme=1).

See rain.f90 to determine the parameter to be defined in callphys.def as a function of the precipitation scheme used.

Physiq.F90 now written in a matricial (more F90) way.
Radii (H2O and CO2 cloud particles, aerosols, duts, ...) calculations now centralized in module radii_mod.F90

and work with the new aerosol scheme implemented by Laura K. Some inconsistency may remain in callsedim.

Implementation compiled with ifort and pgf90.
gcm.e runs in Earth and Early Mars case with CO2 and H2O cycle + dust.

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

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

Download in other formats:

Original Format