Context Navigation

aeropacity.F @ 476

Last change on this file since 476 was 420, checked in by tnavarro, 14 years ago

24/11/11 == TN

corrected minor bug in updatereffrad.F : reffdust was not saved

ccn_factor as not correctly used in sedimentation.

It is now initialized in inifis.F, declared in tracer.h and
used in both simpleclouds.F & callsedim.F to update ice radius.

Commented diagfi outputs in aeropacity.F & improvedclouds.F for non scavenging users.

File size: 17.9 KB

Line
1	SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
2	& pq,tauscaling,tauref,tau,aerosol,reffrad,nueffrad,
3	& QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d)
4
5	! to use 'getin'
6	USE ioipsl_getincom
7	IMPLICIT NONE
8	c=======================================================================
9	c subject:
10	c --------
11	c Computing aerosol optical depth in each gridbox.
12	c
13	c author: F.Forget
14	c ------
15	c update F. Montmessin (water ice scheme)
16	c and S. Lebonnois (12/06/2003) compatibility dust/ice/chemistry
17	c update J.-B. Madeleine 2008-2009:
18	c - added 3D scattering by aerosols;
19	c - dustopacity transferred from physiq.F to callradite.F,
20	c and renamed into aeropacity.F;
21	c
22	c input:
23	c -----
24	c ngrid Number of gridpoint of horizontal grid
25	c nlayer Number of layer
26	c nq Number of tracer
27	c zday Date (time since Ls=0, in martian days)
28	c ls Solar longitude (Ls) , radian
29	c pplay,pplev pressure (Pa) in the middle and boundary of each layer
30	c pq Dust mixing ratio (used if tracer =T and active=T).
31	c reffrad(ngrid,nlayer,naerkind) Aerosol effective radius
32	c QREFvis3d(ngridmx,nlayermx,naerkind) \ 3d extinction coefficients
33	c QREFir3d(ngridmx,nlayermx,naerkind) / at reference wavelengths;
34	c omegaREFvis3d(ngridmx,nlayermx,naerkind) \ 3d single scat. albedo
35	c omegaREFir3d(ngridmx,nlayermx,naerkind) / at reference wavelengths;
36	c
37	c output:
38	c -------
39	c tauref Prescribed mean column optical depth at 700 Pa
40	c tau Column total visible dust optical depth at each point
41	c aerosol aerosol(ig,l,1) is the dust optical
42	c depth in layer l, grid point ig
43
44	c
45	c=======================================================================
46	#include "dimensions.h"
47	#include "dimphys.h"
48	#include "callkeys.h"
49	#include "comcstfi.h"
50	#include "comgeomfi.h"
51	#include "dimradmars.h"
52	#include "yomaer.h"
53	#include "tracer.h"
54	#include "planete.h"
55	#include "aerkind.h"
56
57	c-----------------------------------------------------------------------
58	c
59	c Declarations :
60	c --------------
61	c
62	c Input/Output
63	c ------------
64	INTEGER ngrid,nlayer,nq
65
66	REAL ls,zday,expfactor
67	REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
68	REAL pq(ngrid,nlayer,nq)
69	REAL tauref(ngrid), tau(ngrid,naerkind)
70	REAL aerosol(ngrid,nlayer,naerkind)
71	REAL dsodust(ngridmx,nlayermx)
72	REAL reffrad(ngrid,nlayer,naerkind)
73	REAL nueffrad(ngrid,nlayer,naerkind)
74	REAL QREFvis3d(ngridmx,nlayermx,naerkind)
75	REAL QREFir3d(ngridmx,nlayermx,naerkind)
76	REAL omegaREFvis3d(ngridmx,nlayermx,naerkind)
77	REAL omegaREFir3d(ngridmx,nlayermx,naerkind)
78	c
79	c Local variables :
80	c -----------------
81	INTEGER l,ig,iq,i,j
82	INTEGER iaer ! Aerosol index
83	real topdust(ngridmx)
84	real zlsconst, zp
85	real taueq,tauS,tauN
86	c Mean Qext(vis)/Qext(ir) profile
87	real msolsir(nlayermx,naerkind)
88	c Mean Qext(ir)/Qabs(ir) profile
89	real mqextsqabs(nlayermx,naerkind)
90	c Variables used when multiple particle sizes are used
91	c for dust or water ice particles in the radiative transfer
92	c (see callradite.F for more information).
93	REAL taudusttmp(ngridmx)! Temporary dust opacity
94	! used before scaling
95	REAL tauscaling(ngridmx) ! Scaling factor for qdust and Ndust
96	REAL taudustvis(ngridmx) ! Dust opacity after scaling
97	REAL taudusttes(ngridmx) ! Dust opacity at IR ref. wav. as
98	! "seen" by the GCM.
99	REAL taucloudvis(ngridmx)! Cloud opacity at visible
100	! reference wavelength
101	REAL taucloudtes(ngridmx)! Cloud opacity at infrared
102	! reference wavelength using
103	! Qabs instead of Qext
104	! (direct comparison with TES)
105
106	c local saved variables
107	c ---------------------
108
109	REAL topdust0(ngridmx)
110	SAVE topdust0
111	c Level under which the dust mixing ratio is held constant
112	c when computing the dust opacity in each layer
113	c (this applies when doubleq and active are true)
114	INTEGER, PARAMETER :: cstdustlevel = 7
115
116	LOGICAL firstcall
117	DATA firstcall/.true./
118	SAVE firstcall
119
120	! indexes of water ice and dust tracers:
121	INTEGER,SAVE :: nqdust(nqmx) ! to store the indexes of dust tracers
122	INTEGER,SAVE :: i_ice=0 ! water ice
123	CHARACTER(LEN=20) :: txt ! to temporarly store text
124	CHARACTER(LEN=1) :: txt2 ! to temporarly store text
125	! indexes of dust scatterers:
126	INTEGER,SAVE :: iaerdust(naerkind)
127	INTEGER,SAVE :: naerdust ! number of dust scatterers
128
129	tau(1:ngrid,1:naerkind)=0
130
131	! identify tracers
132
133	IF (firstcall) THEN
134	! identify scatterers that are dust
135	naerdust=0
136	DO iaer=1,naerkind
137	txt=name_iaer(iaer)
138	IF (txt(1:4).eq."dust") THEN
139	naerdust=naerdust+1
140	iaerdust(naerdust)=iaer
141	ENDIF
142	ENDDO
143	! identify tracers which are dust
144	i=0
145	DO iq=1,nq
146	txt=noms(iq)
147	IF (txt(1:4).eq."dust") THEN
148	i=i+1
149	nqdust(i)=iq
150	ENDIF
151	ENDDO
152
153	IF (water.AND.activice) THEN
154	i_ice=igcm_h2o_ice
155	write(,) "aeropacity: i_ice=",i_ice
156	ENDIF
157
158	c altitude of the top of the aerosol layer (km) at Ls=2.76rad:
159	c in the Viking year scenario
160	DO ig=1,ngrid
161	topdust0(ig)=60. -22.SIN(lati(ig))*2
162	END DO
163
164	c typical profile of solsir and (1-w)^(-1):
165	msolsir(1:nlayer,1:naerkind)=0
166	mqextsqabs(1:nlayer,1:naerkind)=0
167	WRITE(,) "Typical profiles of Qext(vis)/Qext(IR)"
168	WRITE(,) " and Qext(IR)/Qabs(IR):"
169	DO iaer = 1, naerkind ! Loop on aerosol kind
170	WRITE(,) "Aerosol # ",iaer
171	DO l=1,nlayer
172	DO ig=1,ngridmx
173	msolsir(l,iaer)=msolsir(l,iaer)+
174	& QREFvis3d(ig,l,iaer)/
175	& QREFir3d(ig,l,iaer)
176	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)+
177	& (1.E0-omegaREFir3d(ig,l,iaer))**(-1)
178	ENDDO
179	msolsir(l,iaer)=msolsir(l,iaer)/REAL(ngridmx)
180	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)/REAL(ngridmx)
181	ENDDO
182	WRITE(,) "solsir: ",msolsir(:,iaer)
183	WRITE(,) "Qext/Qabs(IR): ",mqextsqabs(:,iaer)
184	ENDDO
185
186	! load value of tauvis from callphys.def (if given there,
187	! otherwise default value read from starfi.nc file will be used)
188	call getin("tauvis",tauvis)
189
190	firstcall=.false.
191
192	END IF
193
194	c Vertical column optical depth at 700.Pa
195	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
196	IF(iaervar.eq.1) THEN
197	do ig=1, ngridmx
198	tauref(ig)=max(tauvis,1.e-9) ! tauvis=cste (set in callphys.def
199	! or read in starfi
200	end do
201	ELSE IF (iaervar.eq.2) THEN ! << "Viking" Scenario>>
202
203	tauref(1) = 0.7+.3cos(ls+80.pi/180.) ! like seen by VL1
204	do ig=2,ngrid
205	tauref(ig) = tauref(1)
206	end do
207
208	ELSE IF (iaervar.eq.3) THEN ! << "MGS" scenario >>
209
210	taueq= 0.2 +(0.5-0.2) (cos(0.5(ls-4.363)))**14
211	tauS= 0.1 +(0.5-0.1) (cos(0.5(ls-4.363)))**14
212	tauN = 0.1
213	c if (peri_day.eq.150) then
214	c tauS=0.1
215	c tauN=0.1 +(0.5-0.1) (cos(0.5(ls+pi-4.363)))**14
216	c taueq= 0.2 +(0.5-0.2) (cos(0.5(ls+pi-4.363)))**14
217	c endif
218	do ig=1,ngrid/2 ! Northern hemisphere
219	tauref(ig)= tauN +
220	& (taueq-tauN)0.5(1+tanh((45-lati(ig)180./pi)6/60))
221	end do
222	do ig=ngrid/2+1, ngridmx ! Southern hemisphere
223	tauref(ig)= tauS +
224	& (taueq-tauS)0.5(1+tanh((45+lati(ig)180./pi)6/60))
225	end do
226	ELSE IF ((iaervar.eq.4).or.
227	& ((iaervar.ge.24).and.(iaervar.le.26)))
228	& THEN ! << "TES assimilated dust scenarios >>
229	call readtesassim(ngrid,nlayer,zday,pplev,tauref)
230
231	ELSE IF (iaervar.eq.5) THEN ! << Escalier Scenario>>
232	c tauref(1) = 0.2
233	c if ((ls.ge.210.pi/180.).and.(ls.le.330.pi/180.))
234	c & tauref(1) = 2.5
235	tauref(1) = 2.5
236	if ((ls.ge.30.pi/180.).and.(ls.le.150.pi/180.))
237	& tauref(1) = .2
238
239	do ig=2,ngrid
240	tauref(ig) = tauref(1)
241	end do
242	ELSE
243	stop 'problem with iaervar in aeropacity.F'
244	ENDIF
245
246	c -----------------------------------------------------------------
247	c Computing the opacity in each layer
248	c -----------------------------------------------------------------
249
250	DO iaer = 1, naerkind ! Loop on aerosol kind
251	c --------------------------------------------
252	aerkind: SELECT CASE (name_iaer(iaer))
253	c==================================================================
254	CASE("dust_conrath") aerkind ! Typical dust profile
255	c==================================================================
256
257	c Altitude of the top of the dust layer
258	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
259	zlsconst=SIN(ls-2.76)
260	if (iddist.eq.1) then
261	do ig=1,ngrid
262	topdust(ig)=topdustref ! constant dust layer top
263	end do
264
265	else if (iddist.eq.2) then ! "Viking" scenario
266	do ig=1,ngrid
267	topdust(ig)=topdust0(ig)+18.*zlsconst
268	end do
269
270	else if(iddist.eq.3) then !"MGS" scenario
271	do ig=1,ngrid
272	topdust(ig)=60.+18.*zlsconst
273	& -(32+18zlsconst)sin(lati(ig))**4
274	& - 8zlsconst(sin(lati(ig)))**5
275	end do
276	endif
277
278	c Optical depth in each layer :
279	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
280	if(iddist.ge.1) then
281
282	expfactor=0.
283	DO l=1,nlayer
284	DO ig=1,ngrid
285	c Typical mixing ratio profile
286	if(pplay(ig,l).gt.700.
287	$ /(988.**(topdust(ig)/70.))) then
288	zp=(700./pplay(ig,l))**(70./topdust(ig))
289	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
290	else
291	expfactor=1.e-3
292	endif
293	c Vertical scaling function
294	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) *
295	& expfactor *
296	& QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer)
297	ENDDO
298	ENDDO
299
300	else if(iddist.eq.0) then
301	c old dust vertical distribution function (pollack90)
302	DO l=1,nlayer
303	DO ig=1,ngrid
304	zp=700./pplay(ig,l)
305	aerosol(ig,l,1)= tauref(ig)/700. *
306	s (pplev(ig,l)-pplev(ig,l+1))
307	s max( exp(.03(1.-max(zp,1.))) , 1.E-3 )
308	ENDDO
309	ENDDO
310	end if
311
312	c==================================================================
313	CASE("dust_doubleq") aerkind! Two-moment scheme for dust
314	c (transport of mass and number mixing ratio)
315	c==================================================================
316
317	DO l=1,nlayer
318	IF (l.LE.cstdustlevel) THEN
319	c Opacity in the first levels is held constant to
320	c avoid unrealistic values due to constant lifting:
321	DO ig=1,ngrid
322	aerosol(ig,l,iaer) =
323	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
324	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
325	& pq(ig,cstdustlevel,igcm_dust_mass) *
326	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
327	ENDDO
328	ELSE
329	DO ig=1,ngrid
330	aerosol(ig,l,iaer) =
331	& ( 0.75 * QREFvis3d(ig,l,iaer) /
332	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
333	& pq(ig,l,igcm_dust_mass) *
334	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
335	ENDDO
336	ENDIF
337	ENDDO
338
339	c==================================================================
340	CASE("dust_submicron") aerkind ! Small dust population
341	c==================================================================
342
343	DO l=1,nlayer
344	IF (l.LE.cstdustlevel) THEN
345	c Opacity in the first levels is held constant to
346	c avoid unrealistic values due to constant lifting:
347	DO ig=1,ngrid
348	aerosol(ig,l,iaer) =
349	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
350	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
351	& pq(ig,cstdustlevel,igcm_dust_submicron) *
352	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
353	ENDDO
354	ELSE
355	DO ig=1,ngrid
356	aerosol(ig,l,iaer) =
357	& ( 0.75 * QREFvis3d(ig,l,iaer) /
358	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
359	& pq(ig,l,igcm_dust_submicron) *
360	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
361	ENDDO
362	ENDIF
363	ENDDO
364
365	c==================================================================
366	CASE("h2o_ice") aerkind ! Water ice crystals
367	c==================================================================
368
369	c 1. Initialization
370	aerosol(1:ngrid,1:nlayer,iaer) = 0.
371	taucloudvis(1:ngrid) = 0.
372	taucloudtes(1:ngrid) = 0.
373	c 2. Opacity calculation
374	DO ig=1, ngrid
375	DO l=1,nlayer
376	aerosol(ig,l,iaer) = max(1E-20,
377	& ( 0.75 * QREFvis3d(ig,l,iaer) /
378	& ( rho_ice * reffrad(ig,l,iaer) ) ) *
379	& pq(ig,l,i_ice) *
380	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
381	& )
382	taucloudvis(ig) = taucloudvis(ig) + aerosol(ig,l,iaer)
383	taucloudtes(ig) = taucloudtes(ig) + aerosol(ig,l,iaer)*
384	& QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer) *
385	& ( 1.E0 - omegaREFir3d(ig,l,iaer) )
386	ENDDO
387	ENDDO
388	c 3. Outputs
389	IF (ngrid.NE.1) THEN
390	CALL WRITEDIAGFI(ngridmx,'tauVIS','tauext VIS refwvl',
391	& ' ',2,taucloudvis)
392	CALL WRITEDIAGFI(ngridmx,'tauTES','tauabs IR refwvl',
393	& ' ',2,taucloudtes)
394	IF (callstats) THEN
395	CALL wstats(ngridmx,'tauVIS','tauext VIS refwvl',
396	& ' ',2,taucloudvis)
397	CALL wstats(ngridmx,'tauTES','tauabs IR refwvl',
398	& ' ',2,taucloudtes)
399	ENDIF
400	ELSE
401	c CALL writeg1d(ngrid,1,taucloudtes,'tautes','NU')
402	ENDIF
403	c==================================================================
404	END SELECT aerkind
405	c -----------------------------------
406	ENDDO ! iaer (loop on aerosol kind)
407
408	c -----------------------------------------------------------------
409	c Rescaling each layer to reproduce the choosen (or assimilated)
410	c dust extinction opacity at visible reference wavelength, which
411	c is originally scaled to an equivalent 700Pa pressure surface.
412	c -----------------------------------------------------------------
413
414	c Temporary scaling factor
415	taudusttmp(1:ngrid)=0.
416	DO iaer=1,naerdust
417	DO l=1,nlayer
418	DO ig=1,ngrid
419	c Scaling factor
420	taudusttmp(ig) = taudusttmp(ig) +
421	& aerosol(ig,l,iaerdust(iaer))
422	ENDDO
423	ENDDO
424	ENDDO
425
426	c Saved scaling factor
427	DO ig=1,ngrid
428	tauscaling(ig) = tauref(ig) *
429	& pplev(ig,1) / 700.E0 / taudusttmp(ig)
430	c tauscaling(ig) = 1.e-4
431	ENDDO
432
433	c Opacity computation
434	DO iaer=1,naerdust
435	DO l=1,nlayer
436	DO ig=1,ngrid
437	aerosol(ig,l,iaerdust(iaer)) = max(1E-20,
438	& aerosol(ig,l,iaerdust(iaer))* tauscaling(ig))
439	ENDDO
440	ENDDO
441	ENDDO
442
443	c output for debug
444	c IF (ngrid.NE.1) THEN
445	c CALL WRITEDIAGFI(ngridmx,'taudusttmp','virtual tau dust',
446	c & '#',2,taudusttmp)
447	c CALL WRITEDIAGFI(ngridmx,'tausca','tauscaling',
448	c & '#',2,tauscaling)
449	c ELSE
450	c CALL WRITEDIAGFI(ngridmx,'taudusttmp','virtual tau dust',
451	c & '#',0,taudusttmp)
452	c CALL WRITEDIAGFI(ngridmx,'tausca','tauscaling',
453	c & '#',0,tauscaling)
454	c ENDIF
455	c -----------------------------------------------------------------
456	c Column integrated visible optical depth in each point
457	c -----------------------------------------------------------------
458	DO iaer=1,naerkind
459	do l=1,nlayer
460	do ig=1,ngrid
461	tau(ig,iaer) = tau(ig,iaer) + aerosol(ig,l,iaer)
462	end do
463	end do
464	ENDDO
465	c -----------------------------------------------------------------
466	c Density scaled opacity and column opacity output
467	c -----------------------------------------------------------------
468	c dsodust(1:ngrid,1:nlayer) = 0.
469	c DO iaer=1,naerdust
470	c DO l=1,nlayermx
471	c DO ig=1,ngrid
472	c dsodust(ig,l) = dsodust(ig,l) +
473	c & aerosol(ig,l,iaerdust(iaer)) * g /
474	c & (pplev(ig,l) - pplev(ig,l+1))
475	c ENDDO
476	c ENDDO
477	c IF (ngrid.NE.1) THEN
478	c write(txt2,'(i1.1)') iaer
479	c call WRITEDIAGFI(ngridmx,'taudust'//txt2,
480	c & 'Dust col opacity',
481	c & ' ',2,tau(1,iaerdust(iaer)))
482	c IF (callstats) THEN
483	c CALL wstats(ngridmx,'taudust'//txt2,
484	c & 'Dust col opacity',
485	c & ' ',2,tau(1,iaerdust(iaer)))
486	c ENDIF
487	c ENDIF
488	c ENDDO
489
490	c IF (ngrid.NE.1) THEN
491	c CALL WRITEDIAGFI(ngridmx,'dsodust','tau*g/dp',
492	c & 'm2.kg-1',3,dsodust)
493	c IF (callstats) THEN
494	c CALL wstats(ngridmx,'dsodust',
495	c & 'tau*g/dp',
496	c & 'm2.kg-1',3,dsodust)
497	c ENDIF
498	c ELSE
499	c CALL WRITEDIAGFI(ngrid,"dsodust","dsodust","m2.kg-1",1,
500	c & dsodust)
501	c ENDIF ! of IF (ngrid.NE.1)
502	c -----------------------------------------------------------------
503	return
504	end

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

Download in other formats:

Original Format