Context Navigation

aeropacity.F @ 630

Last change on this file since 630 was 607, checked in by emillour, 13 years ago

Mars GCM:

Updated model so that 1) reference pressure for opacity is now 610 Pa (and
not 700 Pa anymore) and 2) the new MY24-MY30 dust scenarios (input files
dust_MY24.nc, dust_MY25.nc, ..., dust_MY30.nc) can be used
(when setting iaervar=24,25,...,30 in callphys.def): changed
"readtesassim.F90" into "read_dust_scenario.F90" and apadpted aeropacity.F.
one can still use the "old" MY24-MY25-MY26 (input files dust_tes_MY24.nc,
dust_tes_MY25.nc and dust_tes_MY26.nc) scenarios by setting
iaervar=124, 125 or 126 in callphys.def.

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

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

Download in other formats:

Original Format