Context Navigation

aeropacity.F @ 80

Last change on this file since 80 was 38, checked in by emillour, 14 years ago
Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs. EM
File size: 19.3 KB

Line
1	SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
2	& pq,ccn,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 taudustvis(ngridmx) ! Dust opacity after scaling
96	REAL taudusttes(ngridmx) ! Dust opacity at IR ref. wav. as
97	! "seen" by the GCM.
98	REAL taucloudvis(ngridmx)! Cloud opacity at visible
99	! reference wavelength
100	REAL taucloudtes(ngridmx)! Cloud opacity at infrared
101	! reference wavelength using
102	! Qabs instead of Qext
103	! (direct comparison with TES)
104	REAL qdust(ngridmx,nlayermx) ! True dust mass mixing ratio
105	REAL ccn(ngridmx,nlayermx) ! Cloud condensation nuclei
106	! (particules kg-1)
107	REAL qtot(ngridmx) ! Dust column (kg m-2)
108	c
109	c local saved variables
110	c ---------------------
111
112	REAL topdust0(ngridmx)
113	SAVE topdust0
114	c Level under which the dust mixing ratio is held constant
115	c when computing the dust opacity in each layer
116	c (this applies when doubleq and active are true)
117	INTEGER, PARAMETER :: cstdustlevel = 7
118
119	LOGICAL firstcall
120	DATA firstcall/.true./
121	SAVE firstcall
122
123	! indexes of water ice and dust tracers:
124	INTEGER,SAVE :: nqdust(nqmx) ! to store the indexes of dust tracers
125	INTEGER,SAVE :: i_ice=0 ! water ice
126	CHARACTER(LEN=20) :: txt ! to temporarly store text
127	CHARACTER(LEN=1) :: txt2 ! to temporarly store text
128	! indexes of dust scatterers:
129	INTEGER,SAVE :: iaerdust(naerkind)
130	INTEGER,SAVE :: naerdust ! number of dust scatterers
131
132	tau(1:ngrid,1:naerkind)=0
133
134	! identify tracers
135
136	IF (firstcall) THEN
137	! identify scatterers that are dust
138	naerdust=0
139	DO iaer=1,naerkind
140	txt=name_iaer(iaer)
141	IF (txt(1:4).eq."dust") THEN
142	naerdust=naerdust+1
143	iaerdust(naerdust)=iaer
144	ENDIF
145	ENDDO
146	! identify tracers which are dust
147	i=0
148	DO iq=1,nq
149	txt=noms(iq)
150	IF (txt(1:4).eq."dust") THEN
151	i=i+1
152	nqdust(i)=iq
153	ENDIF
154	ENDDO
155
156	IF (water.AND.activice) THEN
157	i_ice=igcm_h2o_ice
158	write(,) "aeropacity: i_ice=",i_ice
159	ENDIF
160
161	c altitude of the top of the aerosol layer (km) at Ls=2.76rad:
162	c in the Viking year scenario
163	DO ig=1,ngrid
164	topdust0(ig)=60. -22.SIN(lati(ig))*2
165	END DO
166
167	c typical profile of solsir and (1-w)^(-1):
168	msolsir(1:nlayer,1:naerkind)=0
169	mqextsqabs(1:nlayer,1:naerkind)=0
170	WRITE(,) "Typical profiles of solsir and Qext/Qabs(IR):"
171	DO iaer = 1, naerkind ! Loop on aerosol kind
172	WRITE(,) "Aerosol # ",iaer
173	DO l=1,nlayer
174	DO ig=1,ngridmx
175	msolsir(l,iaer)=msolsir(l,iaer)+
176	& QREFvis3d(ig,l,iaer)/
177	& QREFir3d(ig,l,iaer)
178	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)+
179	& (1.E0-omegaREFir3d(ig,l,iaer))**(-1)
180	ENDDO
181	msolsir(l,iaer)=msolsir(l,iaer)/REAL(ngridmx)
182	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)/REAL(ngridmx)
183	ENDDO
184	WRITE(,) "solsir: ",msolsir(:,iaer)
185	WRITE(,) "Qext/Qabs(IR): ",mqextsqabs(:,iaer)
186	ENDDO
187
188	! load value of tauvis from callphys.def (if given there,
189	! otherwise default value read from starfi.nc file will be used)
190	call getin("tauvis",tauvis)
191
192	firstcall=.false.
193
194	END IF
195
196	c Vertical column optical depth at 700.Pa
197	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
198	IF(iaervar.eq.1) THEN
199	do ig=1, ngridmx
200	tauref(ig)=max(tauvis,1.e-9) ! tauvis=cste (set in callphys.def
201	! or read in starfi
202	end do
203	ELSE IF (iaervar.eq.2) THEN ! << "Viking" Scenario>>
204
205	tauref(1) = 0.7+.3cos(ls+80.pi/180.) ! like seen by VL1
206	do ig=2,ngrid
207	tauref(ig) = tauref(1)
208	end do
209
210	ELSE IF (iaervar.eq.3) THEN ! << "MGS" scenario >>
211
212	taueq= 0.2 +(0.5-0.2) (cos(0.5(ls-4.363)))**14
213	tauS= 0.1 +(0.5-0.1) (cos(0.5(ls-4.363)))**14
214	tauN = 0.1
215	c if (peri_day.eq.150) then
216	c tauS=0.1
217	c tauN=0.1 +(0.5-0.1) (cos(0.5(ls+pi-4.363)))**14
218	c taueq= 0.2 +(0.5-0.2) (cos(0.5(ls+pi-4.363)))**14
219	c endif
220	do ig=1,ngrid/2 ! Northern hemisphere
221	tauref(ig)= tauN +
222	& (taueq-tauN)0.5(1+tanh((45-lati(ig)180./pi)6/60))
223	end do
224	do ig=ngrid/2+1, ngridmx ! Southern hemisphere
225	tauref(ig)= tauS +
226	& (taueq-tauS)0.5(1+tanh((45+lati(ig)180./pi)6/60))
227	end do
228	ELSE IF ((iaervar.eq.4).or.
229	& ((iaervar.ge.24).and.(iaervar.le.26)))
230	& THEN ! << "TES assimilated dust scenarios >>
231	call readtesassim(ngrid,nlayer,zday,pplev,tauref)
232
233	ELSE IF (iaervar.eq.5) THEN ! << Escalier Scenario>>
234	c tauref(1) = 0.2
235	c if ((ls.ge.210.pi/180.).and.(ls.le.330.pi/180.))
236	c & tauref(1) = 2.5
237	tauref(1) = 2.5
238	if ((ls.ge.30.pi/180.).and.(ls.le.150.pi/180.))
239	& tauref(1) = .2
240
241	do ig=2,ngrid
242	tauref(ig) = tauref(1)
243	end do
244	ELSE
245	stop 'problem with iaervar in aeropacity.F'
246	ENDIF
247
248	c -----------------------------------------------------------------
249	c Computing the opacity in each layer
250	c -----------------------------------------------------------------
251
252	DO iaer = 1, naerkind ! Loop on aerosol kind
253	c --------------------------------------------
254	aerkind: SELECT CASE (name_iaer(iaer))
255	c==================================================================
256	CASE("dust_conrath") aerkind ! Typical dust profile
257	c==================================================================
258
259	c Altitude of the top of the dust layer
260	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
261	zlsconst=SIN(ls-2.76)
262	if (iddist.eq.1) then
263	do ig=1,ngrid
264	topdust(ig)=topdustref ! constant dust layer top
265	end do
266
267	else if (iddist.eq.2) then ! "Viking" scenario
268	do ig=1,ngrid
269	topdust(ig)=topdust0(ig)+18.*zlsconst
270	end do
271
272	else if(iddist.eq.3) then !"MGS" scenario
273	do ig=1,ngrid
274	topdust(ig)=60.+18.*zlsconst
275	& -(32+18zlsconst)sin(lati(ig))**4
276	& - 8zlsconst(sin(lati(ig)))**5
277	end do
278	endif
279
280	c Optical depth in each layer :
281	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
282	if(iddist.ge.1) then
283
284	expfactor=0.
285	DO l=1,nlayer
286	DO ig=1,ngrid
287	c Typical mixing ratio profile
288	if(pplay(ig,l).gt.700.
289	$ /(988.**(topdust(ig)/70.))) then
290	zp=(700./pplay(ig,l))**(70./topdust(ig))
291	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
292	else
293	expfactor=1.e-3
294	endif
295	c Vertical scaling function
296	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) *
297	& expfactor *
298	& QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer)
299	ENDDO
300	ENDDO
301
302	else if(iddist.eq.0) then
303	c old dust vertical distribution function (pollack90)
304	DO l=1,nlayer
305	DO ig=1,ngrid
306	zp=700./pplay(ig,l)
307	aerosol(ig,l,1)= tauref(ig)/700. *
308	s (pplev(ig,l)-pplev(ig,l+1))
309	s max( exp(.03(1.-max(zp,1.))) , 1.E-3 )
310	ENDDO
311	ENDDO
312	end if
313
314	c==================================================================
315	CASE("dust_doubleq") aerkind! Two-moment scheme for dust
316	c (transport of mass and number mixing ratio)
317	c==================================================================
318
319	DO l=1,nlayer
320	IF (l.LE.cstdustlevel) THEN
321	c Opacity in the first levels is held constant to
322	c avoid unrealistic values due to constant lifting:
323	DO ig=1,ngrid
324	aerosol(ig,l,iaer) =
325	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
326	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
327	& pq(ig,cstdustlevel,igcm_dust_mass) *
328	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
329	ENDDO
330	ELSE
331	DO ig=1,ngrid
332	aerosol(ig,l,iaer) =
333	& ( 0.75 * QREFvis3d(ig,l,iaer) /
334	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
335	& pq(ig,l,igcm_dust_mass) *
336	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
337	ENDDO
338	ENDIF
339	ENDDO
340
341	c==================================================================
342	CASE("dust_submicron") aerkind ! Small dust population
343	c==================================================================
344
345	DO l=1,nlayer
346	IF (l.LE.cstdustlevel) THEN
347	c Opacity in the first levels is held constant to
348	c avoid unrealistic values due to constant lifting:
349	DO ig=1,ngrid
350	aerosol(ig,l,iaer) =
351	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
352	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
353	& pq(ig,cstdustlevel,igcm_dust_submicron) *
354	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
355	ENDDO
356	ELSE
357	DO ig=1,ngrid
358	aerosol(ig,l,iaer) =
359	& ( 0.75 * QREFvis3d(ig,l,iaer) /
360	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
361	& pq(ig,l,igcm_dust_submicron) *
362	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
363	ENDDO
364	ENDIF
365	ENDDO
366
367	c==================================================================
368	CASE("h2o_ice") aerkind ! Water ice crystals
369	c==================================================================
370
371	c 1. Initialization
372	aerosol(1:ngrid,1:nlayer,iaer) = 0.
373	taucloudvis(1:ngrid) = 0.
374	taucloudtes(1:ngrid) = 0.
375	c 2. Opacity calculation
376	DO ig=1, ngrid
377	DO l=1,nlayer
378	aerosol(ig,l,iaer) = max(1E-20,
379	& ( 0.75 * QREFvis3d(ig,l,iaer) /
380	& ( rho_ice * reffrad(ig,l,iaer) ) ) *
381	& pq(ig,l,i_ice) *
382	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
383	& )
384	taucloudvis(ig) = taucloudvis(ig) + aerosol(ig,l,iaer)
385	taucloudtes(ig) = taucloudtes(ig) + aerosol(ig,l,iaer)*
386	& QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer) *
387	& ( 1.E0 - omegaREFir3d(ig,l,iaer) )
388	ENDDO
389	ENDDO
390	c 3. Outputs
391	IF (ngrid.NE.1) THEN
392	CALL WRITEDIAGFI(ngridmx,'tauVIS','tauext VIS refwvl',
393	& ' ',2,taucloudvis)
394	CALL WRITEDIAGFI(ngridmx,'tauTES','tauabs IR refwvl',
395	& ' ',2,taucloudtes)
396	IF (callstats) THEN
397	CALL wstats(ngridmx,'tauVIS','tauext VIS refwvl',
398	& ' ',2,taucloudvis)
399	CALL wstats(ngridmx,'tauTES','tauabs IR refwvl',
400	& ' ',2,taucloudtes)
401	ENDIF
402	ELSE
403	c CALL writeg1d(ngrid,1,taucloudtes,'tautes','NU')
404	ENDIF
405	c==================================================================
406	END SELECT aerkind
407	c -----------------------------------
408	ENDDO ! iaer (loop on aerosol kind)
409
410	c -----------------------------------------------------------------
411	c Rescaling each layer to reproduce the choosen (or assimilated)
412	c dust extinction opacity at visible reference wavelength, which
413	c is originally scaled to an equivalent 700Pa pressure surface.
414	c -----------------------------------------------------------------
415
416	taudusttmp(1:ngrid)=0.
417	DO iaer=1,naerdust
418	DO l=1,nlayer
419	DO ig=1,ngrid
420	c Scaling factor
421	taudusttmp(ig) = taudusttmp(ig) +
422	& aerosol(ig,l,iaerdust(iaer))
423	ENDDO
424	ENDDO
425	ENDDO
426	DO iaer=1,naerdust
427	DO l=1,nlayer
428	DO ig=1,ngrid
429	aerosol(ig,l,iaerdust(iaer)) = max(1E-20,
430	& tauref(ig) *
431	& pplev(ig,1) / 700.E0 *
432	& aerosol(ig,l,iaerdust(iaer)) /
433	& taudusttmp(ig)
434	& )
435	ENDDO
436	ENDDO
437	ENDDO
438
439	c -----------------------------------------------------------------
440	c Computing the number of condensation nuclei
441	c -----------------------------------------------------------------
442	DO iaer = 1, naerkind ! Loop on aerosol kind
443	c --------------------------------------------
444	aerkind2: SELECT CASE (name_iaer(iaer))
445	c==================================================================
446	CASE("dust_conrath") aerkind2 ! Typical dust profile
447	c==================================================================
448	DO l=1,nlayer
449	DO ig=1,ngrid
450	ccn(ig,l) = max(aerosol(ig,l,iaer) /
451	& pi / QREFvis3d(ig,l,iaer) *
452	& (1.+nueffrad(ig,l,iaer))**3. /
453	& reffrad(ig,l,iaer)*2. g /
454	& (pplev(ig,l)-pplev(ig,l+1)),1e-30)
455	ENDDO
456	ENDDO
457	c==================================================================
458	CASE("dust_doubleq") aerkind2!Two-moment scheme for dust
459	c (transport of mass and number mixing ratio)
460	c==================================================================
461	qtot(1:ngridmx) = 0.
462	DO l=1,nlayer
463	DO ig=1,ngrid
464	qdust(ig,l) = pq(ig,l,igcm_dust_mass) * tauref(ig) *
465	& pplev(ig,1) / 700.E0 / taudusttmp(ig)
466	qtot(ig) = qtot(ig) + qdust(ig,l) *
467	& (pplev(ig,l)-pplev(ig,l+1)) / g
468	ccn(ig,l) = max( ( ref_r0 /
469	& reffrad(ig,l,iaer) )*3.
470	& r3n_q * qdust(ig,l) ,1e-30)
471	ENDDO
472	ENDDO
473	c==================================================================
474	END SELECT aerkind2
475	c -----------------------------------
476	ENDDO ! iaer (loop on aerosol kind)
477
478	c -----------------------------------------------------------------
479	c Column integrated visible optical depth in each point
480	c -----------------------------------------------------------------
481	DO iaer=1,naerkind
482	do l=1,nlayer
483	do ig=1,ngrid
484	tau(ig,iaer) = tau(ig,iaer) + aerosol(ig,l,iaer)
485	end do
486	end do
487	ENDDO
488	c -----------------------------------------------------------------
489	c Density scaled opacity and column opacity output
490	c -----------------------------------------------------------------
491	dsodust(1:ngrid,1:nlayer) = 0.
492	DO iaer=1,naerdust
493	DO l=1,nlayermx
494	DO ig=1,ngrid
495	dsodust(ig,l) = dsodust(ig,l) +
496	& aerosol(ig,l,iaerdust(iaer)) * g /
497	& (pplev(ig,l) - pplev(ig,l+1))
498	ENDDO
499	ENDDO
500	IF (ngrid.NE.1) THEN
501	write(txt2,'(i1.1)') iaer
502	call WRITEDIAGFI(ngridmx,'taudust'//txt2,
503	& 'Dust col opacity',
504	& ' ',2,tau(1,iaerdust(iaer)))
505	IF (callstats) THEN
506	CALL wstats(ngridmx,'taudust'//txt2,
507	& 'Dust col opacity',
508	& ' ',2,tau(1,iaerdust(iaer)))
509	ENDIF
510	ENDIF
511	ENDDO
512
513	IF (ngrid.NE.1) THEN
514	c CALL WRITEDIAGFI(ngridmx,'dsodust','tau*g/dp',
515	c & 'm2.kg-1',3,dsodust)
516	IF (callstats) THEN
517	CALL wstats(ngridmx,'dsodust',
518	& 'tau*g/dp',
519	& 'm2.kg-1',3,dsodust)
520	ENDIF
521	c CALL WRITEDIAGFI(ngridmx,'ccn','Cond. nuclei',
522	c & 'part kg-1',3,ccn)
523	ELSE
524	CALL writeg1d(ngrid,nlayer,dsodust,'dsodust','m2.kg-1')
525	ENDIF
526	c -----------------------------------------------------------------
527	return
528	end

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

Download in other formats:

Original Format