Context Navigation

aeropacity.F @ 174

Last change on this file since 174 was 83, checked in by aslmd, 14 years ago

mars + LMD_MM_MARS : modifications pour cycle de l'eau, valeurs tunees JBM

used settings reached by JBM to obtain his PhD results

alb_surfice = 0.45 --- in physiq.F and meso_physiq.F
ccn_factor = 4.5 --- in watercloud.F
nuice_sed = 0.45 --- in callsedim.F

NB: this is supposed to be further refined in the future

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

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

Download in other formats:

Original Format