Context Navigation

aeropacity.F @ 390

Last change on this file since 390 was 358, checked in by aslmd, 14 years ago

LMDZ.MARS:

A FIRST VERSION WITH SCAVENGING
TRANSPARENT TO CASUAL USER
FOR WATER CYCLE FOLKS, SEE BELOW

[TN and AS solved conflicts due to r330 >> r357]

07/11/11 == JBM

Changed watercloud.F to call two separate routines,

simpleclouds.F or improvedclouds.F, which are a simplified and
full microphysical scheme for cloud formation, respectively.
Removed the tag called "improved" in watercloud.F, and added
another tag called "microphys" which is defined in callphys.F
instead. Changed routines: callkeys, inifis, physiq, watercloud.

Reimplemented the use of typical profiles of dust particle sizes

and CCNs in simpleclouds.F. Corrected the previously used
analytical CCN profile. Moved ccn_factor to simpleclouds.F,
which won't be used in the new microphysical scheme. Changed
routines: aeropacity, physiq, simpleclouds, watercloud.

Computed rdust at the end of callsedim instead of updatereffrad,

except at firstcall. Renamed rfall into rsedcloud and moved it
in simpleclouds. Moved nuice_sed to initracer.F and added it to
tracer.h. Changed routines: callsedim, physiq, tracer.h,
watercloud, initracer, simpleclouds, updatereffrad.

Added two tracers for the CCNs. Added the different tests in

initracer.F to make sure that, for example, microphys is not
used without doubleq, etc. Corrected an inconsistency in
callsedim.F, and changed the way r0 is updated. Changes
routines: callsedim, inifis, initracer, physiq, testphys1d,
tracer.h.

Added the ability to have a spatially variable density in

newsedim (same method as that used for the radius of
sedimentation). Required the addition of one input to newsedim,
which is the size of the density variable "rho". Changed
routines: callsedim, newsedim.

Added an output to aeropacity called "tauscaling", which is a

factor to convert dust_mass and dust_number to true quantities,
based on tauvis. Removed ccn and qdust from aeropacity, which
are obsolete now.

Wrote improvedcloud.F which includes all the microphysical

scheme, and connected it to the sedimentation scheme. Added and
changed routines: callsedim, physiq, growthrate, nuclea,
improvedclouds, initracer, watercloud, microphys.h.

07/11/11 == TN

Added CCN & dust tracers updates in physiq.F

Corrected a bug that can give negative CCN numbers, due to the
use of single precision values (only 6 decimals) whereas up to 10E+10
CCN can be formed or disappears...

Corrected physical bug that causes h2o_vap or h2o_ice

to be negative in improvedclouds.F.

Corrected physical bug that causes CCN & dust tracers

to be negative in improvedclouds.F when all ice sublimates and
releases dust

Added parameter contact mteta in callphys.def

Default value is still 0.95, see inifis.F

Changed tendancy computation for dust_number in improvedclouds.F

that was not the right one. Indeed, the scavenged dust_number tracer
is computed from the dust_mass one, and its tendancy before scavenging
must be taken into account to compute its scavenging's tendancy.

File size: 17.5 KB

Line
1	SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
2	& pq,tauscaling,tauref,tau,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
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	#include "aerkind.h"
56
57	c-----------------------------------------------------------------------
58	c
59	c Declarations :
60	c --------------
61	c
62	c Input/Output
63	c ------------
64	INTEGER ngrid,nlayer,nq
65
66	REAL ls,zday,expfactor
67	REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
68	REAL pq(ngrid,nlayer,nq)
69	REAL tauref(ngrid), tau(ngrid,naerkind)
70	REAL aerosol(ngrid,nlayer,naerkind)
71	REAL dsodust(ngridmx,nlayermx)
72	REAL reffrad(ngrid,nlayer,naerkind)
73	REAL nueffrad(ngrid,nlayer,naerkind)
74	REAL QREFvis3d(ngridmx,nlayermx,naerkind)
75	REAL QREFir3d(ngridmx,nlayermx,naerkind)
76	REAL omegaREFvis3d(ngridmx,nlayermx,naerkind)
77	REAL omegaREFir3d(ngridmx,nlayermx,naerkind)
78	c
79	c Local variables :
80	c -----------------
81	INTEGER l,ig,iq,i,j
82	INTEGER iaer ! Aerosol index
83	real topdust(ngridmx)
84	real zlsconst, zp
85	real taueq,tauS,tauN
86	c Mean Qext(vis)/Qext(ir) profile
87	real msolsir(nlayermx,naerkind)
88	c Mean Qext(ir)/Qabs(ir) profile
89	real mqextsqabs(nlayermx,naerkind)
90	c Variables used when multiple particle sizes are used
91	c for dust or water ice particles in the radiative transfer
92	c (see callradite.F for more information).
93	REAL taudusttmp(ngridmx)! Temporary dust opacity
94	! used before scaling
95	REAL tauscaling(ngridmx) ! Scaling factor for qdust and Ndust
96	REAL taudustvis(ngridmx) ! Dust opacity after scaling
97	REAL taudusttes(ngridmx) ! Dust opacity at IR ref. wav. as
98	! "seen" by the GCM.
99	REAL taucloudvis(ngridmx)! Cloud opacity at visible
100	! reference wavelength
101	REAL taucloudtes(ngridmx)! Cloud opacity at infrared
102	! reference wavelength using
103	! Qabs instead of Qext
104	! (direct comparison with TES)
105
106	c local saved variables
107	c ---------------------
108
109	REAL topdust0(ngridmx)
110	SAVE topdust0
111	c Level under which the dust mixing ratio is held constant
112	c when computing the dust opacity in each layer
113	c (this applies when doubleq and active are true)
114	INTEGER, PARAMETER :: cstdustlevel = 7
115
116	LOGICAL firstcall
117	DATA firstcall/.true./
118	SAVE firstcall
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	CHARACTER(LEN=20) :: txt ! to temporarly store text
124	CHARACTER(LEN=1) :: txt2 ! to temporarly store text
125	! indexes of dust scatterers:
126	INTEGER,SAVE :: iaerdust(naerkind)
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 altitude of the top of the aerosol layer (km) at Ls=2.76rad:
159	c in the Viking year scenario
160	DO ig=1,ngrid
161	topdust0(ig)=60. -22.SIN(lati(ig))*2
162	END DO
163
164	c typical profile of solsir and (1-w)^(-1):
165	msolsir(1:nlayer,1:naerkind)=0
166	mqextsqabs(1:nlayer,1:naerkind)=0
167	WRITE(,) "Typical profiles of Qext(vis)/Qext(IR)"
168	WRITE(,) " and Qext(IR)/Qabs(IR):"
169	DO iaer = 1, naerkind ! Loop on aerosol kind
170	WRITE(,) "Aerosol # ",iaer
171	DO l=1,nlayer
172	DO ig=1,ngridmx
173	msolsir(l,iaer)=msolsir(l,iaer)+
174	& QREFvis3d(ig,l,iaer)/
175	& QREFir3d(ig,l,iaer)
176	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)+
177	& (1.E0-omegaREFir3d(ig,l,iaer))**(-1)
178	ENDDO
179	msolsir(l,iaer)=msolsir(l,iaer)/REAL(ngridmx)
180	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)/REAL(ngridmx)
181	ENDDO
182	WRITE(,) "solsir: ",msolsir(:,iaer)
183	WRITE(,) "Qext/Qabs(IR): ",mqextsqabs(:,iaer)
184	ENDDO
185
186	! load value of tauvis from callphys.def (if given there,
187	! otherwise default value read from starfi.nc file will be used)
188	call getin("tauvis",tauvis)
189
190	firstcall=.false.
191
192	END IF
193
194	c Vertical column optical depth at 700.Pa
195	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
196	IF(iaervar.eq.1) THEN
197	do ig=1, ngridmx
198	tauref(ig)=max(tauvis,1.e-9) ! tauvis=cste (set in callphys.def
199	! or read in starfi
200	end do
201	ELSE IF (iaervar.eq.2) THEN ! << "Viking" Scenario>>
202
203	tauref(1) = 0.7+.3cos(ls+80.pi/180.) ! like seen by VL1
204	do ig=2,ngrid
205	tauref(ig) = tauref(1)
206	end do
207
208	ELSE IF (iaervar.eq.3) THEN ! << "MGS" scenario >>
209
210	taueq= 0.2 +(0.5-0.2) (cos(0.5(ls-4.363)))**14
211	tauS= 0.1 +(0.5-0.1) (cos(0.5(ls-4.363)))**14
212	tauN = 0.1
213	c if (peri_day.eq.150) then
214	c tauS=0.1
215	c tauN=0.1 +(0.5-0.1) (cos(0.5(ls+pi-4.363)))**14
216	c taueq= 0.2 +(0.5-0.2) (cos(0.5(ls+pi-4.363)))**14
217	c endif
218	do ig=1,ngrid/2 ! Northern hemisphere
219	tauref(ig)= tauN +
220	& (taueq-tauN)0.5(1+tanh((45-lati(ig)180./pi)6/60))
221	end do
222	do ig=ngrid/2+1, ngridmx ! Southern hemisphere
223	tauref(ig)= tauS +
224	& (taueq-tauS)0.5(1+tanh((45+lati(ig)180./pi)6/60))
225	end do
226	ELSE IF ((iaervar.eq.4).or.
227	& ((iaervar.ge.24).and.(iaervar.le.26)))
228	& THEN ! << "TES assimilated dust scenarios >>
229	call readtesassim(ngrid,nlayer,zday,pplev,tauref)
230
231	ELSE IF (iaervar.eq.5) THEN ! << Escalier Scenario>>
232	c tauref(1) = 0.2
233	c if ((ls.ge.210.pi/180.).and.(ls.le.330.pi/180.))
234	c & tauref(1) = 2.5
235	tauref(1) = 2.5
236	if ((ls.ge.30.pi/180.).and.(ls.le.150.pi/180.))
237	& tauref(1) = .2
238
239	do ig=2,ngrid
240	tauref(ig) = tauref(1)
241	end do
242	ELSE
243	stop 'problem with iaervar in aeropacity.F'
244	ENDIF
245
246	c -----------------------------------------------------------------
247	c Computing the opacity in each layer
248	c -----------------------------------------------------------------
249
250	DO iaer = 1, naerkind ! Loop on aerosol kind
251	c --------------------------------------------
252	aerkind: SELECT CASE (name_iaer(iaer))
253	c==================================================================
254	CASE("dust_conrath") aerkind ! Typical dust profile
255	c==================================================================
256
257	c Altitude of the top of the dust layer
258	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
259	zlsconst=SIN(ls-2.76)
260	if (iddist.eq.1) then
261	do ig=1,ngrid
262	topdust(ig)=topdustref ! constant dust layer top
263	end do
264
265	else if (iddist.eq.2) then ! "Viking" scenario
266	do ig=1,ngrid
267	topdust(ig)=topdust0(ig)+18.*zlsconst
268	end do
269
270	else if(iddist.eq.3) then !"MGS" scenario
271	do ig=1,ngrid
272	topdust(ig)=60.+18.*zlsconst
273	& -(32+18zlsconst)sin(lati(ig))**4
274	& - 8zlsconst(sin(lati(ig)))**5
275	end do
276	endif
277
278	c Optical depth in each layer :
279	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
280	if(iddist.ge.1) then
281
282	expfactor=0.
283	DO l=1,nlayer
284	DO ig=1,ngrid
285	c Typical mixing ratio profile
286	if(pplay(ig,l).gt.700.
287	$ /(988.**(topdust(ig)/70.))) then
288	zp=(700./pplay(ig,l))**(70./topdust(ig))
289	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
290	else
291	expfactor=1.e-3
292	endif
293	c Vertical scaling function
294	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) *
295	& expfactor *
296	& QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer)
297	ENDDO
298	ENDDO
299
300	else if(iddist.eq.0) then
301	c old dust vertical distribution function (pollack90)
302	DO l=1,nlayer
303	DO ig=1,ngrid
304	zp=700./pplay(ig,l)
305	aerosol(ig,l,1)= tauref(ig)/700. *
306	s (pplev(ig,l)-pplev(ig,l+1))
307	s max( exp(.03(1.-max(zp,1.))) , 1.E-3 )
308	ENDDO
309	ENDDO
310	end if
311
312	c==================================================================
313	CASE("dust_doubleq") aerkind! Two-moment scheme for dust
314	c (transport of mass and number mixing ratio)
315	c==================================================================
316
317	DO l=1,nlayer
318	IF (l.LE.cstdustlevel) THEN
319	c Opacity in the first levels is held constant to
320	c avoid unrealistic values due to constant lifting:
321	DO ig=1,ngrid
322	aerosol(ig,l,iaer) =
323	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
324	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
325	& pq(ig,cstdustlevel,igcm_dust_mass) *
326	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
327	ENDDO
328	ELSE
329	DO ig=1,ngrid
330	aerosol(ig,l,iaer) =
331	& ( 0.75 * QREFvis3d(ig,l,iaer) /
332	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
333	& pq(ig,l,igcm_dust_mass) *
334	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
335	ENDDO
336	ENDIF
337	ENDDO
338
339	c==================================================================
340	CASE("dust_submicron") aerkind ! Small dust population
341	c==================================================================
342
343	DO l=1,nlayer
344	IF (l.LE.cstdustlevel) THEN
345	c Opacity in the first levels is held constant to
346	c avoid unrealistic values due to constant lifting:
347	DO ig=1,ngrid
348	aerosol(ig,l,iaer) =
349	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
350	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
351	& pq(ig,cstdustlevel,igcm_dust_submicron) *
352	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
353	ENDDO
354	ELSE
355	DO ig=1,ngrid
356	aerosol(ig,l,iaer) =
357	& ( 0.75 * QREFvis3d(ig,l,iaer) /
358	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
359	& pq(ig,l,igcm_dust_submicron) *
360	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
361	ENDDO
362	ENDIF
363	ENDDO
364
365	c==================================================================
366	CASE("h2o_ice") aerkind ! Water ice crystals
367	c==================================================================
368
369	c 1. Initialization
370	aerosol(1:ngrid,1:nlayer,iaer) = 0.
371	taucloudvis(1:ngrid) = 0.
372	taucloudtes(1:ngrid) = 0.
373	c 2. Opacity calculation
374	DO ig=1, ngrid
375	DO l=1,nlayer
376	aerosol(ig,l,iaer) = max(1E-20,
377	& ( 0.75 * QREFvis3d(ig,l,iaer) /
378	& ( rho_ice * reffrad(ig,l,iaer) ) ) *
379	& pq(ig,l,i_ice) *
380	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
381	& )
382	taucloudvis(ig) = taucloudvis(ig) + aerosol(ig,l,iaer)
383	taucloudtes(ig) = taucloudtes(ig) + aerosol(ig,l,iaer)*
384	& QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer) *
385	& ( 1.E0 - omegaREFir3d(ig,l,iaer) )
386	ENDDO
387	ENDDO
388	c 3. Outputs
389	IF (ngrid.NE.1) THEN
390	CALL WRITEDIAGFI(ngridmx,'tauVIS','tauext VIS refwvl',
391	& ' ',2,taucloudvis)
392	CALL WRITEDIAGFI(ngridmx,'tauTES','tauabs IR refwvl',
393	& ' ',2,taucloudtes)
394	IF (callstats) THEN
395	CALL wstats(ngridmx,'tauVIS','tauext VIS refwvl',
396	& ' ',2,taucloudvis)
397	CALL wstats(ngridmx,'tauTES','tauabs IR refwvl',
398	& ' ',2,taucloudtes)
399	ENDIF
400	ELSE
401	c CALL writeg1d(ngrid,1,taucloudtes,'tautes','NU')
402	ENDIF
403	c==================================================================
404	END SELECT aerkind
405	c -----------------------------------
406	ENDDO ! iaer (loop on aerosol kind)
407
408	c -----------------------------------------------------------------
409	c Rescaling each layer to reproduce the choosen (or assimilated)
410	c dust extinction opacity at visible reference wavelength, which
411	c is originally scaled to an equivalent 700Pa pressure surface.
412	c -----------------------------------------------------------------
413
414	c Temporary scaling factor
415	taudusttmp(1:ngrid)=0.
416	DO iaer=1,naerdust
417	DO l=1,nlayer
418	DO ig=1,ngrid
419	c Scaling factor
420	taudusttmp(ig) = taudusttmp(ig) +
421	& aerosol(ig,l,iaerdust(iaer))
422	ENDDO
423	ENDDO
424	ENDDO
425
426	c Saved scaling factor
427	DO ig=1,ngrid
428	tauscaling(ig) = tauref(ig) *
429	& pplev(ig,1) / 700.E0 / taudusttmp(ig)
430	ENDDO
431
432	c Opacity computation
433	DO iaer=1,naerdust
434	DO l=1,nlayer
435	DO ig=1,ngrid
436	aerosol(ig,l,iaerdust(iaer)) = max(1E-20,
437	& tauref(ig) *
438	& pplev(ig,1) / 700.E0 *
439	& aerosol(ig,l,iaerdust(iaer)) /
440	& taudusttmp(ig)
441	& )
442	ENDDO
443	ENDDO
444	ENDDO
445
446	c -----------------------------------------------------------------
447	c Column integrated visible optical depth in each point
448	c -----------------------------------------------------------------
449	DO iaer=1,naerkind
450	do l=1,nlayer
451	do ig=1,ngrid
452	tau(ig,iaer) = tau(ig,iaer) + aerosol(ig,l,iaer)
453	end do
454	end do
455	ENDDO
456	c -----------------------------------------------------------------
457	c Density scaled opacity and column opacity output
458	c -----------------------------------------------------------------
459	dsodust(1:ngrid,1:nlayer) = 0.
460	DO iaer=1,naerdust
461	DO l=1,nlayermx
462	DO ig=1,ngrid
463	dsodust(ig,l) = dsodust(ig,l) +
464	& aerosol(ig,l,iaerdust(iaer)) * g /
465	& (pplev(ig,l) - pplev(ig,l+1))
466	ENDDO
467	ENDDO
468	IF (ngrid.NE.1) THEN
469	write(txt2,'(i1.1)') iaer
470	call WRITEDIAGFI(ngridmx,'taudust'//txt2,
471	& 'Dust col opacity',
472	& ' ',2,tau(1,iaerdust(iaer)))
473	IF (callstats) THEN
474	CALL wstats(ngridmx,'taudust'//txt2,
475	& 'Dust col opacity',
476	& ' ',2,tau(1,iaerdust(iaer)))
477	ENDIF
478	ENDIF
479	ENDDO
480
481	IF (ngrid.NE.1) THEN
482	c CALL WRITEDIAGFI(ngridmx,'dsodust','tau*g/dp',
483	c & 'm2.kg-1',3,dsodust)
484	IF (callstats) THEN
485	CALL wstats(ngridmx,'dsodust',
486	& 'tau*g/dp',
487	& 'm2.kg-1',3,dsodust)
488	ENDIF
489	ELSE
490	CALL WRITEDIAGFI(ngrid,"dsodust","dsodust","m2.kg-1",1,
491	& dsodust)
492	ENDIF ! of IF (ngrid.NE.1)
493	c -----------------------------------------------------------------
494	return
495	end

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

Download in other formats:

Original Format