Context Navigation

aeropacity.F @ 706

Last change on this file since 706 was 677, checked in by emillour, 13 years ago

Mars GCM:

Added the possibility to use the "cold" (iaervar=6) or "warm" (iaervar=7) dust scenarios (derived from the 7 MY24 to MY30 scenarios) in aeropacity.F and read_dust_scenario.F

File size: 18.5 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.6).and.(iaervar.le.7)) THEN
239	! cold or warm synthetic scenarios
240	call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
241	ELSE IF ((iaervar.ge.24).and.(iaervar.le.30))
242	& THEN ! << MY... dust scenarios >>
243	call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
244	ELSE IF ((iaervar.eq.4).or.
245	& ((iaervar.ge.124).and.(iaervar.le.126))) THEN
246	! "old" TES assimation dust scenario (values at 700Pa in files!)
247	call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
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.odpref
293	$ /(988.**(topdust(ig)/70.))) then
294	zp=(odpref/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=odpref/pplay(ig,l)
311	aerosol(ig,l,1)= tauref(ig)/odpref *
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 -- Now done in physiq.F
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	! 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 odpref Pa pressure surface.
418	c -----------------------------------------------------------------
419
420	c Temporary scaling factor
421	taudusttmp(1:ngrid)=0.
422	DO iaer=1,naerdust
423	DO l=1,nlayer
424	DO ig=1,ngrid
425	c Scaling factor
426	taudusttmp(ig) = taudusttmp(ig) +
427	& aerosol(ig,l,iaerdust(iaer))
428	ENDDO
429	ENDDO
430	ENDDO
431
432	c Saved scaling factor
433	DO ig=1,ngrid
434	tauscaling(ig) = tauref(ig) *
435	& pplev(ig,1) / odpref / taudusttmp(ig)
436	c tauscaling(ig) = 1.e-4
437	ENDDO
438
439	c Opacity computation
440	DO iaer=1,naerdust
441	DO l=1,nlayer
442	DO ig=1,ngrid
443	aerosol(ig,l,iaerdust(iaer)) = max(1E-20,
444	& aerosol(ig,l,iaerdust(iaer))* tauscaling(ig))
445	ENDDO
446	ENDDO
447	ENDDO
448
449	c output for debug
450	c IF (ngrid.NE.1) THEN
451	c CALL WRITEDIAGFI(ngridmx,'taudusttmp','virtual tau dust',
452	c & '#',2,taudusttmp)
453	c CALL WRITEDIAGFI(ngridmx,'tausca','tauscaling',
454	c & '#',2,tauscaling)
455	c ELSE
456	c CALL WRITEDIAGFI(ngridmx,'taudusttmp','virtual tau dust',
457	c & '#',0,taudusttmp)
458	c CALL WRITEDIAGFI(ngridmx,'tausca','tauscaling',
459	c & '#',0,tauscaling)
460	c ENDIF
461	c -----------------------------------------------------------------
462	c Column integrated visible optical depth in each point
463	c -----------------------------------------------------------------
464	DO iaer=1,naerkind
465	do l=1,nlayer
466	do ig=1,ngrid
467	tau(ig,iaer) = tau(ig,iaer) + aerosol(ig,l,iaer)
468	end do
469	end do
470	ENDDO
471	c -----------------------------------------------------------------
472	c Density scaled opacity and column opacity output
473	c -----------------------------------------------------------------
474	c dsodust(1:ngrid,1:nlayer) = 0.
475	c DO iaer=1,naerdust
476	c DO l=1,nlayermx
477	c DO ig=1,ngrid
478	c dsodust(ig,l) = dsodust(ig,l) +
479	c & aerosol(ig,l,iaerdust(iaer)) * g /
480	c & (pplev(ig,l) - pplev(ig,l+1))
481	c ENDDO
482	c ENDDO
483	c IF (ngrid.NE.1) THEN
484	c write(txt2,'(i1.1)') iaer
485	c call WRITEDIAGFI(ngridmx,'taudust'//txt2,
486	c & 'Dust col opacity',
487	c & ' ',2,tau(1,iaerdust(iaer)))
488	c IF (callstats) THEN
489	c CALL wstats(ngridmx,'taudust'//txt2,
490	c & 'Dust col opacity',
491	c & ' ',2,tau(1,iaerdust(iaer)))
492	c ENDIF
493	c ENDIF
494	c ENDDO
495
496	c IF (ngrid.NE.1) THEN
497	c CALL WRITEDIAGFI(ngridmx,'dsodust','tau*g/dp',
498	c & 'm2.kg-1',3,dsodust)
499	c IF (callstats) THEN
500	c CALL wstats(ngridmx,'dsodust',
501	c & 'tau*g/dp',
502	c & 'm2.kg-1',3,dsodust)
503	c ENDIF
504	c ELSE
505	c CALL WRITEDIAGFI(ngrid,"dsodust","dsodust","m2.kg-1",1,
506	c & dsodust)
507	c ENDIF ! of IF (ngrid.NE.1)
508	c -----------------------------------------------------------------
509	return
510	end

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

Download in other formats:

Original Format