Context Navigation

aeropacity.F @ 1104

Last change on this file since 1104 was 1104, checked in by tnavarro, 11 years ago
bug correction of the bug correction
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	use comgeomfi_h, only: lati ! grid point latitudes (rad)
10	use yomaer_h, only: tauvis
11	IMPLICIT NONE
12	c=======================================================================
13	c subject:
14	c --------
15	c Computing aerosol optical depth in each gridbox.
16	c
17	c author: F.Forget
18	c ------
19	c update F. Montmessin (water ice scheme)
20	c and S. Lebonnois (12/06/2003) compatibility dust/ice/chemistry
21	c update J.-B. Madeleine 2008-2009:
22	c - added 3D scattering by aerosols;
23	c - dustopacity transferred from physiq.F to callradite.F,
24	c and renamed into aeropacity.F;
25	c update E. Millour, march 2012:
26	c - reference pressure is now set to 610Pa (not 700Pa)
27	c
28	c input:
29	c -----
30	c ngrid Number of gridpoint of horizontal grid
31	c nlayer Number of layer
32	c nq Number of tracer
33	c zday Date (time since Ls=0, in martian days)
34	c ls Solar longitude (Ls) , radian
35	c pplay,pplev pressure (Pa) in the middle and boundary of each layer
36	c pq Dust mixing ratio (used if tracer =T and active=T).
37	c reffrad(ngrid,nlayer,naerkind) Aerosol effective radius
38	c QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients
39	c QREFir3d(ngrid,nlayer,naerkind) / at reference wavelengths;
40	c omegaREFvis3d(ngrid,nlayer,naerkind) \ 3d single scat. albedo
41	c omegaREFir3d(ngrid,nlayer,naerkind) / at reference wavelengths;
42	c
43	c output:
44	c -------
45	c tauref Prescribed mean column optical depth at 610 Pa
46	c tau Column total visible dust optical depth at each point
47	c aerosol aerosol(ig,l,1) is the dust optical
48	c depth in layer l, grid point ig
49
50	c
51	c=======================================================================
52	!#include "dimensions.h"
53	!#include "dimphys.h"
54	#include "callkeys.h"
55	#include "comcstfi.h"
56	!#include "comgeomfi.h"
57	!#include "dimradmars.h"
58	!#include "yomaer.h"
59	!#include "tracer.h"
60	! naerkind is set in scatterers.h (built when compiling with makegcm -s #)
61	#include"scatterers.h"
62	#include "planete.h"
63	#include "aerkind.h"
64
65	c-----------------------------------------------------------------------
66	c
67	c Declarations :
68	c --------------
69	c
70	c Input/Output
71	c ------------
72	INTEGER ngrid,nlayer,nq
73
74	REAL ls,zday,expfactor
75	REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
76	REAL pq(ngrid,nlayer,nq)
77	REAL tauref(ngrid), tau(ngrid,naerkind)
78	REAL aerosol(ngrid,nlayer,naerkind)
79	REAL dsodust(ngrid,nlayer)
80	REAL reffrad(ngrid,nlayer,naerkind)
81	REAL nueffrad(ngrid,nlayer,naerkind)
82	REAL QREFvis3d(ngrid,nlayer,naerkind)
83	REAL QREFir3d(ngrid,nlayer,naerkind)
84	REAL omegaREFvis3d(ngrid,nlayer,naerkind)
85	REAL omegaREFir3d(ngrid,nlayer,naerkind)
86	c
87	c Local variables :
88	c -----------------
89	INTEGER l,ig,iq,i,j
90	INTEGER iaer ! Aerosol index
91	real topdust(ngrid)
92	real zlsconst, zp
93	real taueq,tauS,tauN
94	c Mean Qext(vis)/Qext(ir) profile
95	real msolsir(nlayer,naerkind)
96	c Mean Qext(ir)/Qabs(ir) profile
97	real mqextsqabs(nlayer,naerkind)
98	c Variables used when multiple particle sizes are used
99	c for dust or water ice particles in the radiative transfer
100	c (see callradite.F for more information).
101	REAL taudusttmp(ngrid)! Temporary dust opacity
102	! used before scaling
103	REAL tauscaling(ngrid) ! Scaling factor for qdust and Ndust
104	REAL taudustvis(ngrid) ! Dust opacity after scaling
105	REAL taudusttes(ngrid) ! Dust opacity at IR ref. wav. as
106	! "seen" by the GCM.
107	REAL taucloudvis(ngrid)! Cloud opacity at visible
108	! reference wavelength
109	REAL taucloudtes(ngrid)! Cloud opacity at infrared
110	! reference wavelength using
111	! Qabs instead of Qext
112	! (direct comparison with TES)
113
114	c local saved variables
115	c ---------------------
116
117	REAL,SAVE,ALLOCATABLE :: topdust0(:)
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 :: cstdustlevel = 7
123
124	LOGICAL,SAVE :: firstcall=.true.
125
126	! indexes of water ice and dust tracers:
127	INTEGER,ALLOCATABLE,SAVE :: nqdust(:) ! to store the indexes of 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	! allocate local saved arrays
142	allocate(nqdust(nq))
143	allocate(topdust0(ngrid))
144
145	! identify scatterers that are dust
146	naerdust=0
147	DO iaer=1,naerkind
148	txt=name_iaer(iaer)
149	IF (txt(1:4).eq."dust") THEN
150	naerdust=naerdust+1
151	iaerdust(naerdust)=iaer
152	ENDIF
153	ENDDO
154	! identify tracers which are dust
155	i=0
156	DO iq=1,nq
157	txt=noms(iq)
158	IF (txt(1:4).eq."dust") THEN
159	i=i+1
160	nqdust(i)=iq
161	ENDIF
162	ENDDO
163
164	IF (water.AND.activice) THEN
165	i_ice=igcm_h2o_ice
166	write(,) "aeropacity: i_ice=",i_ice
167	ENDIF
168
169	c altitude of the top of the aerosol layer (km) at Ls=2.76rad:
170	c in the Viking year scenario
171	DO ig=1,ngrid
172	topdust0(ig)=60. -22.SIN(lati(ig))*2
173	END DO
174
175	c typical profile of solsir and (1-w)^(-1):
176	msolsir(1:nlayer,1:naerkind)=0
177	mqextsqabs(1:nlayer,1:naerkind)=0
178	WRITE(,) "Typical profiles of Qext(vis)/Qext(IR)"
179	WRITE(,) " and Qext(IR)/Qabs(IR):"
180	DO iaer = 1, naerkind ! Loop on aerosol kind
181	WRITE(,) "Aerosol # ",iaer
182	DO l=1,nlayer
183	DO ig=1,ngrid
184	msolsir(l,iaer)=msolsir(l,iaer)+
185	& QREFvis3d(ig,l,iaer)/
186	& QREFir3d(ig,l,iaer)
187	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)+
188	& (1.E0-omegaREFir3d(ig,l,iaer))**(-1)
189	ENDDO
190	msolsir(l,iaer)=msolsir(l,iaer)/REAL(ngrid)
191	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)/REAL(ngrid)
192	ENDDO
193	WRITE(,) "solsir: ",msolsir(:,iaer)
194	WRITE(,) "Qext/Qabs(IR): ",mqextsqabs(:,iaer)
195	ENDDO
196
197	! load value of tauvis from callphys.def (if given there,
198	! otherwise default value read from starfi.nc file will be used)
199	call getin("tauvis",tauvis)
200
201	firstcall=.false.
202
203	END IF
204
205	c Vertical column optical depth at "odpref" Pa
206	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
207	IF(freedust) THEN
208	tauref(:) = 0. ! tauref is computed after, instead of being forced
209
210	ELSE IF(iaervar.eq.1) THEN
211	do ig=1, ngrid
212	tauref(ig)=max(tauvis,1.e-9) ! tauvis=cste (set in callphys.def
213	! or read in starfi
214	end do
215	ELSE IF (iaervar.eq.2) THEN ! << "Viking" Scenario>>
216
217	tauref(1) = 0.7+.3cos(ls+80.pi/180.) ! like seen by VL1
218	do ig=2,ngrid
219	tauref(ig) = tauref(1)
220	end do
221
222	ELSE IF (iaervar.eq.3) THEN ! << "MGS" scenario >>
223
224	taueq= 0.2 +(0.5-0.2) (cos(0.5(ls-4.363)))**14
225	tauS= 0.1 +(0.5-0.1) (cos(0.5(ls-4.363)))**14
226	tauN = 0.1
227	c if (peri_day.eq.150) then
228	c tauS=0.1
229	c tauN=0.1 +(0.5-0.1) (cos(0.5(ls+pi-4.363)))**14
230	c taueq= 0.2 +(0.5-0.2) (cos(0.5(ls+pi-4.363)))**14
231	c endif
232	do ig=1,ngrid
233	if (lati(ig).ge.0) then
234	! Northern hemisphere
235	tauref(ig)= tauN +
236	& (taueq-tauN)0.5(1+tanh((45-lati(ig)180./pi)6/60))
237	else
238	! Southern hemisphere
239	tauref(ig)= tauS +
240	& (taueq-tauS)0.5(1+tanh((45+lati(ig)180./pi)6/60))
241	endif
242	enddo ! of do ig=1,ngrid
243	ELSE IF (iaervar.eq.5) THEN ! << Escalier Scenario>>
244	c tauref(1) = 0.2
245	c if ((ls.ge.210.pi/180.).and.(ls.le.330.pi/180.))
246	c & tauref(1) = 2.5
247	tauref(1) = 2.5
248	if ((ls.ge.30.pi/180.).and.(ls.le.150.pi/180.))
249	& tauref(1) = .2
250
251	do ig=2,ngrid
252	tauref(ig) = tauref(1)
253	end do
254	ELSE IF ((iaervar.ge.6).and.(iaervar.le.7)) THEN
255	! cold or warm synthetic scenarios
256	call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
257	ELSE IF ((iaervar.ge.24).and.(iaervar.le.30))
258	& THEN ! << MY... dust scenarios >>
259	call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
260	ELSE IF ((iaervar.eq.4).or.
261	& ((iaervar.ge.124).and.(iaervar.le.126))) THEN
262	! "old" TES assimation dust scenario (values at 700Pa in files!)
263	call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
264	ELSE
265	stop 'problem with iaervar in aeropacity.F'
266	ENDIF
267
268	c -----------------------------------------------------------------
269	c Computing the opacity in each layer
270	c -----------------------------------------------------------------
271
272	DO iaer = 1, naerkind ! Loop on aerosol kind
273	c --------------------------------------------
274	aerkind: SELECT CASE (name_iaer(iaer))
275	c==================================================================
276	CASE("dust_conrath") aerkind ! Typical dust profile
277	c==================================================================
278
279	c Altitude of the top of the dust layer
280	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
281	zlsconst=SIN(ls-2.76)
282	if (iddist.eq.1) then
283	do ig=1,ngrid
284	topdust(ig)=topdustref ! constant dust layer top
285	end do
286
287	else if (iddist.eq.2) then ! "Viking" scenario
288	do ig=1,ngrid
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	c tauscaling(ig) = 1.e-4
457	ENDDO
458
459	ENDIF
460
461	c Opacity computation
462	DO iaer=1,naerdust
463	DO l=1,nlayer
464	DO ig=1,ngrid
465	aerosol(ig,l,iaerdust(iaer)) = max(1E-20,
466	& aerosol(ig,l,iaerdust(iaer))* tauscaling(ig))
467	ENDDO
468	ENDDO
469	ENDDO
470
471	IF (freedust) THEN
472	! tauref has been initialized to 0 before.
473	DO iaer=1,naerdust
474	DO l=1,nlayer
475	DO ig=1,ngrid
476	tauref(ig) = tauref(ig) +
477	& aerosol(ig,l,iaerdust(iaer))
478	ENDDO
479	ENDDO
480	ENDDO
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