Context Navigation

aeropacity.F @ 1242

Last change on this file since 1242 was 1226, checked in by aslmd, 11 years ago
LMDZ.MARS : Replaced comcstfi and planete includes by modules.
File size: 19.4 KB

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

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

Download in other formats:

Original Format