Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

assim_aeropacity.F @ 3094

Last change on this file since 3094 was 42, checked in by aslmd, 14 years ago
LMD_MM_MARS: tests avec physique completement commune avec le GCM nouvelle physique
File size: 16.1 KB

Line
1	SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,pq,
2	& tauref,tau,aerosol,reffrad,
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
22	c input:
23	c -----
24	c ngrid Number of gridpoint of horizontal grid
25	c nlayer Number of layer
26	c nq Number of tracer
27	c zday Date (time since Ls=0, in martian days)
28	c ls Solar longitude (Ls) , radian
29	c pplay,pplev pressure (Pa) in the middle and boundary of each layer
30	c pq Dust mixing ratio (used if tracer =T and active=T).
31	c reffrad(ngrid,nlayer,naerkind) Aerosol effective radius
32	c QREFvis3d(ngridmx,nlayermx,naerkind) \ 3d extinction coefficients
33	c QREFir3d(ngridmx,nlayermx,naerkind) / at reference wavelengths;
34	c omegaREFvis3d(ngridmx,nlayermx,naerkind) \ 3d single scat. albedo
35	c omegaREFir3d(ngridmx,nlayermx,naerkind) / at reference wavelengths;
36	c
37	c output:
38	c -------
39	c tauref Prescribed mean column optical depth at 700 Pa
40	c tau Column total visible dust optical depth at each point
41	c aerosol aerosol(ig,l,1) is the dust optical
42	c depth in layer l, grid point ig
43
44	c
45	c=======================================================================
46	#include "dimensions.h"
47	#include "dimphys.h"
48	#include "callkeys.h"
49	#include "comcstfi.h"
50	#include "comgeomfi.h"
51	#include "dimradmars.h"
52	#include "yomaer.h"
53	#include "tracer.h"
54	#include "planete.h"
55
56	c-----------------------------------------------------------------------
57	c
58	c Declarations :
59	c --------------
60	c
61	c Input/Output
62	c ------------
63	INTEGER ngrid,nlayer,nq
64
65	REAL ls,zday,expfactor
66	REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
67	REAL pq(ngrid,nlayer,nq)
68	REAL tauref(ngrid), tau(ngrid,naerkind)
69	REAL aerosol(ngrid,nlayer,naerkind)
70	REAL reffrad(ngrid,nlayer,naerkind)
71	REAL QREFvis3d(ngridmx,nlayermx,naerkind)
72	REAL QREFir3d(ngridmx,nlayermx,naerkind)
73	REAL omegaREFvis3d(ngridmx,nlayermx,naerkind)
74	REAL omegaREFir3d(ngridmx,nlayermx,naerkind)
75	c
76	c Local variables :
77	c -----------------
78	INTEGER l,ig,iq,i,j
79	INTEGER iaer ! Aerosol index
80	real topdust(ngridmx)
81	real zlsconst, zp
82	real taueq,tauS,tauN
83	real r0,reff,coefsize
84	c Mean Qext(vis)/Qext(ir) profile
85	real msolsir(nlayermx,naerkind)
86	c Mean Qext(ir)/Qabs(ir) profile
87	real mqextsqabs(nlayermx,naerkind)
88	c Variables used when multiple particle sizes are used
89	c for dust or water ice particles in the radiative transfer
90	c (see callradite.F for more information).
91	REAL taudusttmp(ngridmx)! Temporary dust opacity
92	! used before scaling
93	REAL taudustvis(ngridmx) ! Dust opacity after scaling
94	REAL taudusttes(ngridmx) ! Dust opacity at IR ref. wav. as
95	! "seen" by the GCM.
96	REAL taucloudvis(ngridmx)! Cloud opacity at visible
97	! reference wavelength
98	REAL taucloudtes(ngridmx)! Cloud opacity at infrared
99	! reference wavelength using
100	! Qabs instead of Qext
101	! (direct comparison with TES)
102	c
103	c local saved variables
104	c ---------------------
105
106	REAL topdust0(ngridmx)
107	SAVE topdust0
108
109	LOGICAL firstcall
110	DATA firstcall/.true./
111	SAVE firstcall
112
113	! indexes of water ice and dust tracers:
114	INTEGER,SAVE :: nqdust(nqmx) ! to store the indexes of dust tracers
115	INTEGER,SAVE :: i_ice=0 ! water ice
116	CHARACTER(LEN=20) :: tracername ! to temporarly store text
117
118	call zerophys(ngrid*naerkind,tau)
119
120	! identify tracers
121
122	IF (firstcall) THEN
123	! identify tracers which are dust
124	i=0
125	DO iq=1,nq
126	tracername=noms(iq)
127	IF (tracername(1:4).eq."dust") THEN
128	i=i+1
129	nqdust(i)=iq
130	ENDIF
131	ENDDO
132
133	IF (water.AND.activice) THEN
134	i_ice=igcm_h2o_ice
135	write(,) "aeropacity: i_ice=",i_ice
136	ENDIF
137
138	c altitude of the top of the aerosol layer (km) at Ls=2.76rad:
139	c in the Viking year scenario
140	DO ig=1,ngrid
141	topdust0(ig)=60. -22.SIN(lati(ig))*2
142	END DO
143
144	c typical profile of solsir and (1-w)^(-1):
145	call zerophys(nlayer*naerkind,msolsir)
146	call zerophys(nlayer*naerkind,mqextsqabs)
147	WRITE(,) "Typical profiles of solsir and Qext/Qabs(IR):"
148	DO iaer = 1, naerkind ! Loop on aerosol kind
149	WRITE(,) "Aerosol # ",iaer
150	DO l=1,nlayer
151	DO ig=1,ngridmx
152	msolsir(l,iaer)=msolsir(l,iaer)+
153	& QREFvis3d(ig,l,iaer)/
154	& QREFir3d(ig,l,iaer)
155	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)+
156	& (1.E0-omegaREFir3d(ig,l,iaer))**(-1)
157	ENDDO
158	msolsir(l,iaer)=msolsir(l,iaer)/REAL(ngridmx)
159	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)/REAL(ngridmx)
160	ENDDO
161	WRITE(,) "solsir: ",msolsir(:,iaer)
162	WRITE(,) "Qext/Qabs(IR): ",mqextsqabs(:,iaer)
163	ENDDO
164
165	! load value of tauvis from callphys.def (if given there,
166	! otherwise default value read from starfi.nc file will be used)
167	call getin("tauvis",tauvis)
168
169	firstcall=.false.
170
171	END IF
172
173	DO iaer = 1, naerkind ! Loop on aerosol kind
174	c --------------------------------------------
175	aerkind: SELECT CASE (iaer)
176	c==================================================================
177	CASE(1) aerkind ! Dust (iaer=1)
178	c==================================================================
179
180	c -------------------------------------------------------------
181	c 1) Prescribed dust (if tracer=F or active=F)
182	c -------------------------------------------------------------
183	IF ((.not.tracer) .or. (.not.active)) THEN
184
185	c Vertical column optical depth at 700.Pa
186	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
187	IF(iaervar.eq.1) THEN
188	do ig=1, ngridmx
189	tauref(ig)=max(tauvis,1.e-9) ! tauvis=cste (set in callphys.def
190	! or read in starfi
191	end do
192	ELSE IF (iaervar.eq.2) THEN ! << "Viking" Scenario>>
193
194	tauref(1) = 0.7+.3cos(ls+80.pi/180.) ! like seen by VL1
195	do ig=2,ngrid
196	tauref(ig) = tauref(1)
197	end do
198
199	ELSE IF (iaervar.eq.3) THEN ! << "MGS" scenario >>
200
201	taueq= 0.2 +(0.5-0.2) (cos(0.5(ls-4.363)))**14
202	tauS= 0.1 +(0.5-0.1) (cos(0.5(ls-4.363)))**14
203	tauN = 0.1
204	c if (peri_day.eq.150) then
205	c tauS=0.1
206	c tauN=0.1 +(0.5-0.1) (cos(0.5(ls+pi-4.363)))**14
207	c taueq= 0.2 +(0.5-0.2) (cos(0.5(ls+pi-4.363)))**14
208	c endif
209	do ig=1,ngrid/2 ! Northern hemisphere
210	tauref(ig)= tauN +
211	& (taueq-tauN)0.5(1+tanh((45-lati(ig)180./pi)6/60))
212	end do
213	do ig=ngrid/2+1, ngridmx ! Southern hemisphere
214	tauref(ig)= tauS +
215	& (taueq-tauS)0.5(1+tanh((45+lati(ig)180./pi)6/60))
216	end do
217
218	ELSE IF ((iaervar.eq.4).or.
219	& ((iaervar.ge.24).and.(iaervar.le.26)))
220	& THEN ! << "TES assimilated dust scenarios >>
221	call readtesassim(ngrid,nlayer,zday,pplev,tauref)
222
223	ELSE IF (iaervar.eq.5) THEN ! << Escalier Scenario>>
224	c tauref(1) = 0.2
225	c if ((ls.ge.210.pi/180.).and.(ls.le.330.pi/180.))
226	c & tauref(1) = 2.5
227	tauref(1) = 2.5
228	if ((ls.ge.30.pi/180.).and.(ls.le.150.pi/180.))
229	& tauref(1) = .2
230	do ig=2,ngrid
231	tauref(ig) = tauref(1)
232	end do
233
234	ELSE IF (iaervar.gt.99) THEN ! << input netcdf file >>
235	c*************WRF
236	c
237	c 2. customized dust opacity field
238	c ex: from assimilation
239	call meso_readtesassim(ngrid,nlayer,zday,pplev,tauref,
240	. iaervar)
241	c
242	c*************WRF
243
244	ELSE
245	stop 'problem with iaervar in aeropacity.F'
246	ENDIF
247
248	c Altitude of the top of the dust layer
249	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
250	zlsconst=SIN(ls-2.76)
251	if (iddist.eq.1) then
252	do ig=1,ngrid
253	topdust(ig)=topdustref ! constant dust layer top
254	end do
255
256	else if (iddist.eq.2) then ! "Viking" scenario
257	do ig=1,ngrid
258	topdust(ig)=topdust0(ig)+18.*zlsconst
259	end do
260
261	else if(iddist.eq.3) then !"MGS" scenario
262	do ig=1,ngrid
263	topdust(ig)=60.+18.*zlsconst
264	& -(32+18zlsconst)sin(lati(ig))**4
265	& - 8zlsconst(sin(lati(ig)))**5
266	end do
267	endif
268
269
270	c Optical depth in each layer :
271	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
272	if(iddist.ge.1) then
273	expfactor=0.
274	CALL zerophys(ngrid,taudusttmp)
275	DO l=1,nlayer
276	DO ig=1,ngrid
277	c Typical mixing ratio profile
278	if(pplay(ig,l).gt.700.
279	$ /(988.**(topdust(ig)/70.))) then
280	zp=(700./pplay(ig,l))**(70./topdust(ig))
281	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
282	else
283	expfactor=1.e-3
284	endif
285	c Vertical scaling function
286	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) *
287	& expfactor *
288	& QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer)
289	c Scaling factor
290	taudusttmp(ig)=taudusttmp(ig)+aerosol(ig,l,iaer)
291	ENDDO
292	ENDDO
293
294	c Rescaling each layer to reproduce the choosen (or
295	c assimilated) dust extinction opacity at visible
296	c reference wavelength, which is originally scaled
297	c to an equivalent 700Pa pressure surface.
298	DO l=1,nlayer
299	DO ig=1,ngrid
300	aerosol(ig,l,iaer) = tauref(ig) *
301	& pplev(ig,1) / 700.E0 *
302	& aerosol(ig,l,iaer) / taudusttmp(ig)
303	ENDDO
304	ENDDO
305
306	CALL zerophys(ngrid,taudustvis)
307	CALL zerophys(ngrid,taudusttes)
308	DO l=1,nlayer
309	DO ig=1,ngrid
310	taudustvis(ig) = taudustvis(ig) + aerosol(ig,l,iaer)
311	taudusttes(ig) = taudusttes(ig) + aerosol(ig,l,iaer)*
312	& QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer)*
313	& ( 1.E0 - omegaREFir3d(ig,l,iaer) )
314	ENDDO
315	ENDDO
316
317	c Outputs
318	IF (ngrid.NE.1) THEN
319	CALL WRITEDIAGFI(ngridmx,'taudustTES','dust abs IR',
320	& ' ',2,taudusttes)
321	CALL wstats(ngridmx,'taudustTES','dust abs IR',
322	& ' ',2,taudusttes)
323	ELSE
324	CALL writeg1d(ngrid,1,taudusttes,'taudusttes','NU')
325	ENDIF
326
327	c changement dans le calcul de la distribution verticale
328	c dans le cas des scenarios de poussieres assimiles
329	c if (iaervar.eq.4) THEN ! TES
330	c call zerophys(ngrid*naerkind,tau)
331	c
332	c do l=1,nlayer
333	c do ig=1,ngrid
334	c tau(ig,1)=tau(ig,1)+ aerosol(ig,l,1)
335	c end do
336	c end do
337	c do l=1,nlayer
338	c do ig=1,ngrid
339	c aerosol(ig,l,1)=aerosol(ig,l,1)*tauref(ig)/tau(ig,1)
340	c $ *(pplev(ig,1)/700)
341	c end do
342	c end do
343	c endif
344	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
345	else if(iddist.eq.0) then
346	c old dust vertical distribution function (pollack90)
347	DO l=1,nlayer
348	DO ig=1,ngrid
349	zp=700./pplay(ig,l)
350	aerosol(ig,l,1)= tauref(ig)/700. *
351	s (pplev(ig,l)-pplev(ig,l+1))
352	s max( exp(.03(1.-max(zp,1.))) , 1.E-3 )
353	ENDDO
354	ENDDO
355	end if
356
357	c ---------------------------------------------------------------------
358	c 2) Transported radiatively active dust (if tracer=T and active=T)
359	c ----------------------------------------------------------------------
360	ELSE IF ((tracer) .and. (active)) THEN
361	c The dust opacity is computed from q
362
363	c a) "doubleq" technique (transport of mass and number mixing ratio)
364	c ~~~~~~~~~~~~~~~~~~~
365	if(doubleq) then
366
367	call zerophys(ngridnlayernaerkind,aerosol)
368
369	c Computing effective radius :
370	do l=1,nlayer
371	do ig=1, ngrid
372	r0=
373	& (r3n_qpq(ig,l,1)/max(pq(ig,l,2),0.01))*(1./3.)
374	r0=min(max(r0,1.e-10),500.e-6)
375	reff= ref_r0 * r0
376	cc If reff is small, the transported dust mean Qext
377	c is reduced from the reference dust Qext by a factor "coefsize"
378
379	coefsize=min(max(2.52e6*reff-0.043 ,0.) ,1.)
380
381	cc It is added 1.e-8 to pq to avoid low
382
383	aerosol(ig,l,1)=aerosol(ig,l,1)+ 1.E-8 +
384	& ( 0.75Qext(1)coefsize/(rho_dust*reff))
385	& * (pq(ig,l,nqdust(1)))*
386	c only one dust bin to use with doubleq
387	& (pplev(ig,l)-pplev(ig,l+1))/g
388	end do
389	end do
390	call zerophys(ngrid,tauref)
391
392	c b) Size bin technique (each aerosol can contribute to opacity))
393	c ~~~~~~~~~~~~~~~~~~
394	else
395	c The dust opacity is computed from q
396	call zerophys(ngridnlayernaerkind,aerosol)
397	do iq=1,dustbin
398	do l=1,nlayer
399	do ig=1,ngrid
400	cc qextrhor(iq) is (3/4)Qext/(rhoreff)
401	cc It is added 1.e-8 to pq to avoid low
402	aerosol(ig,l,1)=aerosol(ig,l,1)+
403	& qextrhor(nqdust(iq))(pq(ig,l,nqdust(iq))+1.e-8)
404	& (pplev(ig,l)-pplev(ig,l+1))/g
405	end do
406	end do
407	end do
408	call zerophys(ngrid,tauref)
409	end if ! (doubleq)
410	END IF ! (dust scenario)
411
412
413	c==================================================================
414	CASE(2) aerkind ! Water ice crystals (iaer=2)
415	c==================================================================
416
417	IF (water.AND.activice) THEN
418	c 1. Initialization
419	CALL zerophys(ngrid*nlayer,aerosol(1,1,iaer))
420	CALL zerophys(ngrid,taucloudvis)
421	CALL zerophys(ngrid,taucloudtes)
422	c 2. Opacity calculation
423	DO ig=1, ngrid
424	DO l=1,nlayer
425	aerosol(ig,l,iaer) =
426	& ( 0.75 * QREFvis3d(ig,l,iaer) /
427	& ( rho_ice * reffrad(ig,l,iaer) ) ) *
428	& ( pq(ig,l,i_ice) + 1.E-8 ) *
429	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
430	taucloudvis(ig) = taucloudvis(ig) + aerosol(ig,l,iaer)
431	taucloudtes(ig) = taucloudtes(ig) + aerosol(ig,l,iaer)*
432	& QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer) *
433	& ( 1.E0 - omegaREFir3d(ig,l,iaer) )
434	ENDDO
435	ENDDO
436	c 3. Outputs
437	IF (ngrid.NE.1) THEN
438	CALL WRITEDIAGFI(ngridmx,'tauTES','tauabs IR refwvl',
439	& ' ',2,taucloudtes)
440	CALL wstats(ngridmx,'tauTES','tauabs IR refwvl',
441	& ' ',2,taucloudtes)
442	ELSE
443	CALL writeg1d(ngrid,1,taucloudtes,'tautes','NU')
444	ENDIF
445	ENDIF
446
447	c==================================================================
448	END SELECT aerkind
449
450	c -----------------------------------------------------------------
451	c Column integrated visible optical depth in each point
452	c -----------------------------------------------------------------
453
454	do l=1,nlayer
455	do ig=1,ngrid
456	tau(ig,iaer) = tau(ig,iaer) + aerosol(ig,l,iaer)
457	end do
458	end do
459
460	c -----------------------------------
461	ENDDO ! iaer (loop on aerosol kind)
462
463	return
464	end

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

Download in other formats:

Original Format