Context Navigation

aeropacity.F @ 1292

Last change on this file since 1292 was 1278, checked in by emillour, 11 years ago

Mars GCM:

Add possibility to use clim or MY31 dust scenarios (and realistic EUV for MY31).

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

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

Download in other formats:

Original Format