Context Navigation

aeropacity.F @ 1502

Last change on this file since 1502 was 1502, checked in by emillour, 9 years ago

Mars GCM:

added missing condition to use MY32 dust scenario in aeropacity
corrected stats output for vmr_ice*rice in physiq

File size: 25.8 KB

Line
1	SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
2	& pq,tauscaling,tauref,tau,taucloudtes,aerosol,dsodust,reffrad,
3	& nueffrad,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	& nqdust
10	use comgeomfi_h, only: lati, sinlat ! grid point latitudes (rad)
11	#ifdef DUSTSTORM
12	use comgeomfi_h, only: long
13	use tracer_mod, only: r3n_q, ref_r0, igcm_dust_number
14	#endif
15	use planete_h
16	USE comcstfi_h
17	use dimradmars_mod, only: naerkind, name_iaer,
18	& iaerdust,tauvis,
19	& iaer_dust_conrath,iaer_dust_doubleq,
20	& iaer_dust_submicron,iaer_h2o_ice
21	IMPLICIT NONE
22	c=======================================================================
23	c subject:
24	c --------
25	c Computing aerosol optical depth in each gridbox.
26	c
27	c author: F.Forget
28	c ------
29	c update F. Montmessin (water ice scheme)
30	c and S. Lebonnois (12/06/2003) compatibility dust/ice/chemistry
31	c update J.-B. Madeleine 2008-2009:
32	c - added 3D scattering by aerosols;
33	c - dustopacity transferred from physiq.F to callradite.F,
34	c and renamed into aeropacity.F;
35	c update E. Millour, march 2012:
36	c - reference pressure is now set to 610Pa (not 700Pa)
37	c
38	c input:
39	c -----
40	c ngrid Number of gridpoint of horizontal grid
41	c nlayer Number of layer
42	c nq Number of tracer
43	c zday Date (time since Ls=0, in martian days)
44	c ls Solar longitude (Ls) , radian
45	c pplay,pplev pressure (Pa) in the middle and boundary of each layer
46	c pq Dust mixing ratio (used if tracer =T and active=T).
47	c reffrad(ngrid,nlayer,naerkind) Aerosol effective radius
48	c QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients
49	c QREFir3d(ngrid,nlayer,naerkind) / at reference wavelengths;
50	c omegaREFvis3d(ngrid,nlayer,naerkind) \ 3d single scat. albedo
51	c omegaREFir3d(ngrid,nlayer,naerkind) / at reference wavelengths;
52	c
53	c output:
54	c -------
55	c tauref Prescribed mean column optical depth at 610 Pa
56	c tau Column total visible dust optical depth at each point
57	c aerosol aerosol(ig,l,1) is the dust optical
58	c depth in layer l, grid point ig
59
60	c
61	c=======================================================================
62	#include "callkeys.h"
63
64	c-----------------------------------------------------------------------
65	c
66	c Declarations :
67	c --------------
68	c
69	c Input/Output
70	c ------------
71	INTEGER ngrid,nlayer,nq
72
73	REAL ls,zday,expfactor
74	REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
75	REAL pq(ngrid,nlayer,nq)
76	REAL tauref(ngrid), tau(ngrid,naerkind)
77	REAL aerosol(ngrid,nlayer,naerkind)
78	REAL dsodust(ngrid,nlayer)
79	REAL reffrad(ngrid,nlayer,naerkind)
80	REAL nueffrad(ngrid,nlayer,naerkind)
81	REAL QREFvis3d(ngrid,nlayer,naerkind)
82	REAL QREFir3d(ngrid,nlayer,naerkind)
83	REAL omegaREFvis3d(ngrid,nlayer,naerkind)
84	REAL omegaREFir3d(ngrid,nlayer,naerkind)
85	c
86	c Local variables :
87	c -----------------
88	INTEGER l,ig,iq,i,j
89	INTEGER iaer ! Aerosol index
90	real topdust(ngrid)
91	real zlsconst, zp
92	real taueq,tauS,tauN
93	c Mean Qext(vis)/Qext(ir) profile
94	real msolsir(nlayer,naerkind)
95	c Mean Qext(ir)/Qabs(ir) profile
96	real mqextsqabs(nlayer,naerkind)
97	c Variables used when multiple particle sizes are used
98	c for dust or water ice particles in the radiative transfer
99	c (see callradite.F for more information).
100	REAL taudusttmp(ngrid)! Temporary dust opacity
101	! used before scaling
102	REAL tauscaling(ngrid) ! Scaling factor for qdust and Ndust
103	REAL taudustvis(ngrid) ! Dust opacity after scaling
104	REAL taudusttes(ngrid) ! Dust opacity at IR ref. wav. as
105	! "seen" by the GCM.
106	REAL taucloudvis(ngrid)! Cloud opacity at visible
107	! reference wavelength
108	REAL taucloudtes(ngrid)! Cloud opacity at infrared
109	! reference wavelength using
110	! Qabs instead of Qext
111	! (direct comparison with TES)
112	REAL topdust0(ngrid)
113
114	#ifdef DUSTSTORM
115	!! Local dust storms
116	logical localstorm ! =true to create a local dust storm
117	real taulocref,ztoploc,radloc,lonloc,latloc ! local dust storm parameters
118	real reffstorm, yeah
119	REAL ray(ngrid) ! distance from dust storm center
120	REAL tauuser(ngrid) ! opacity perturbation due to dust storm
121	REAL more_dust(ngrid,nlayer,2) ! Mass mixing ratio perturbation due to the dust storm
122	REAL int_factor(ngrid) ! useful factor to compute mmr perturbation
123	real l_top ! layer of the storm's top
124	REAL zalt(ngrid, nlayer) ! useful factor to compute l_top
125	#endif
126
127	c local saved variables
128	c ---------------------
129
130	c Level under which the dust mixing ratio is held constant
131	c when computing the dust opacity in each layer
132	c (this applies when doubleq and active are true)
133	INTEGER, PARAMETER :: cstdustlevel0 = 7
134	INTEGER, SAVE :: cstdustlevel
135
136	LOGICAL,SAVE :: firstcall=.true.
137
138	! indexes of water ice and dust tracers:
139	INTEGER,SAVE :: i_ice=0 ! water ice
140	real,parameter :: odpref=610. ! DOD reference pressure (Pa)
141	CHARACTER(LEN=20) :: txt ! to temporarly store text
142	CHARACTER(LEN=1) :: txt2 ! to temporarly store text
143	! indexes of dust scatterers:
144	INTEGER,SAVE :: naerdust ! number of dust scatterers
145
146	tau(1:ngrid,1:naerkind)=0
147
148	! identify tracers
149
150	IF (firstcall) THEN
151	! identify scatterers that are dust
152	naerdust=0
153	DO iaer=1,naerkind
154	txt=name_iaer(iaer)
155	IF (txt(1:4).eq."dust") THEN
156	naerdust=naerdust+1
157	iaerdust(naerdust)=iaer
158	ENDIF
159	ENDDO
160	! identify tracers which are dust
161	i=0
162	DO iq=1,nq
163	txt=noms(iq)
164	IF (txt(1:4).eq."dust") THEN
165	i=i+1
166	nqdust(i)=iq
167	ENDIF
168	ENDDO
169
170	IF (water.AND.activice) THEN
171	i_ice=igcm_h2o_ice
172	write(,) "aeropacity: i_ice=",i_ice
173	ENDIF
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	IF (freedust) THEN
202	cstdustlevel = 1
203	ELSE
204	cstdustlevel = cstdustlevel0
205	ENDIF
206
207
208	#ifndef DUSTSTORM
209	firstcall=.false.
210	#endif
211
212	END IF
213
214	c Vertical column optical depth at "odpref" Pa
215	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
216	IF(freedust) THEN
217	tauref(:) = 0. ! tauref is computed after, instead of being forced
218
219	ELSE IF(iaervar.eq.1) THEN
220	do ig=1, ngrid
221	tauref(ig)=max(tauvis,1.e-9) ! tauvis=cste (set in callphys.def
222	! or read in starfi
223	end do
224	ELSE IF (iaervar.eq.2) THEN ! << "Viking" Scenario>>
225
226	tauref(1) = 0.7+.3cos(ls+80.pi/180.) ! like seen by VL1
227	do ig=2,ngrid
228	tauref(ig) = tauref(1)
229	end do
230
231	ELSE IF (iaervar.eq.3) THEN ! << "MGS" scenario >>
232
233	taueq= 0.2 +(0.5-0.2) (cos(0.5(ls-4.363)))**14
234	tauS= 0.1 +(0.5-0.1) (cos(0.5(ls-4.363)))**14
235	tauN = 0.1
236	c if (peri_day.eq.150) then
237	c tauS=0.1
238	c tauN=0.1 +(0.5-0.1) (cos(0.5(ls+pi-4.363)))**14
239	c taueq= 0.2 +(0.5-0.2) (cos(0.5(ls+pi-4.363)))**14
240	c endif
241	do ig=1,ngrid
242	if (lati(ig).ge.0) then
243	! Northern hemisphere
244	tauref(ig)= tauN +
245	& (taueq-tauN)0.5(1+tanh((45-lati(ig)180./pi)6/60))
246	else
247	! Southern hemisphere
248	tauref(ig)= tauS +
249	& (taueq-tauS)0.5(1+tanh((45+lati(ig)180./pi)6/60))
250	endif
251	enddo ! of do ig=1,ngrid
252	ELSE IF (iaervar.eq.5) THEN ! << Escalier Scenario>>
253	c tauref(1) = 0.2
254	c if ((ls.ge.210.pi/180.).and.(ls.le.330.pi/180.))
255	c & tauref(1) = 2.5
256	tauref(1) = 2.5
257	if ((ls.ge.30.pi/180.).and.(ls.le.150.pi/180.))
258	& tauref(1) = .2
259
260	do ig=2,ngrid
261	tauref(ig) = tauref(1)
262	end do
263	ELSE IF ((iaervar.ge.6).and.(iaervar.le.8)) THEN
264	! clim, cold or warm synthetic scenarios
265	call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
266	ELSE IF ((iaervar.ge.24).and.(iaervar.le.32))
267	& THEN ! << MY... dust scenarios >>
268	call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
269	ELSE IF ((iaervar.eq.4).or.
270	& ((iaervar.ge.124).and.(iaervar.le.126))) THEN
271	! "old" TES assimation dust scenario (values at 700Pa in files!)
272	call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
273	ELSE
274	stop 'problem with iaervar in aeropacity.F'
275	ENDIF
276
277	c -----------------------------------------------------------------
278	c Computing the opacity in each layer
279	c -----------------------------------------------------------------
280
281	DO iaer = 1, naerkind ! Loop on aerosol kind
282	c --------------------------------------------
283	aerkind: SELECT CASE (name_iaer(iaer))
284	c==================================================================
285	CASE("dust_conrath") aerkind ! Typical dust profile
286	c==================================================================
287
288	c Altitude of the top of the dust layer
289	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
290	zlsconst=SIN(ls-2.76)
291	if (iddist.eq.1) then
292	do ig=1,ngrid
293	topdust(ig)=topdustref ! constant dust layer top
294	end do
295
296	else if (iddist.eq.2) then ! "Viking" scenario
297	do ig=1,ngrid
298	! altitude of the top of the aerosol layer (km) at Ls=2.76rad:
299	! in the Viking year scenario
300	topdust0(ig)=60. -22.sinlat(ig)*2
301	topdust(ig)=topdust0(ig)+18.*zlsconst
302	end do
303
304	else if(iddist.eq.3) then !"MGS" scenario
305	do ig=1,ngrid
306	topdust(ig)=60.+18.*zlsconst
307	& -(32+18zlsconst)sin(lati(ig))**4
308	& - 8zlsconst(sin(lati(ig)))**5
309	end do
310	endif
311
312	c Optical depth in each layer :
313	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
314	if(iddist.ge.1) then
315
316	expfactor=0.
317	DO l=1,nlayer
318	DO ig=1,ngrid
319	c Typical mixing ratio profile
320	if(pplay(ig,l).gt.odpref
321	$ /(988.**(topdust(ig)/70.))) then
322	zp=(odpref/pplay(ig,l))**(70./topdust(ig))
323	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
324	else
325	expfactor=1.e-3
326	endif
327	c Vertical scaling function
328	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) *
329	& expfactor *
330	& QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer)
331	ENDDO
332	ENDDO
333
334	else if(iddist.eq.0) then
335	c old dust vertical distribution function (pollack90)
336	DO l=1,nlayer
337	DO ig=1,ngrid
338	zp=odpref/pplay(ig,l)
339	aerosol(ig,l,1)= tauref(ig)/odpref *
340	s (pplev(ig,l)-pplev(ig,l+1))
341	s max( exp(.03(1.-max(zp,1.))) , 1.E-3 )
342	ENDDO
343	ENDDO
344	end if
345
346	c==================================================================
347	CASE("dust_doubleq") aerkind! Two-moment scheme for dust
348	c (transport of mass and number mixing ratio)
349	c==================================================================
350
351	DO l=1,nlayer
352	IF (l.LE.cstdustlevel) THEN
353	c Opacity in the first levels is held constant to
354	c avoid unrealistic values due to constant lifting:
355	DO ig=1,ngrid
356	aerosol(ig,l,iaer) =
357	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
358	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
359	& pq(ig,cstdustlevel,igcm_dust_mass) *
360	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
361	! DENSITY SCALED OPACITY IN INFRARED:
362	dsodust(ig,l) =
363	& ( 0.75 * QREFir3d(ig,cstdustlevel,iaer) /
364	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
365	& pq(ig,cstdustlevel,igcm_dust_mass)
366	ENDDO
367	ELSE
368	DO ig=1,ngrid
369	aerosol(ig,l,iaer) =
370	& ( 0.75 * QREFvis3d(ig,l,iaer) /
371	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
372	& pq(ig,l,igcm_dust_mass) *
373	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
374	! DENSITY SCALED OPACITY IN INFRARED:
375	dsodust(ig,l) =
376	& ( 0.75 * QREFir3d(ig,l,iaer) /
377	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
378	& pq(ig,l,igcm_dust_mass)
379	ENDDO
380	ENDIF
381	ENDDO
382
383	c==================================================================
384	CASE("dust_submicron") aerkind ! Small dust population
385	c==================================================================
386
387	DO l=1,nlayer
388	IF (l.LE.cstdustlevel) THEN
389	c Opacity in the first levels is held constant to
390	c avoid unrealistic values due to constant lifting:
391	DO ig=1,ngrid
392	aerosol(ig,l,iaer) =
393	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
394	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
395	& pq(ig,cstdustlevel,igcm_dust_submicron) *
396	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
397	ENDDO
398	ELSE
399	DO ig=1,ngrid
400	aerosol(ig,l,iaer) =
401	& ( 0.75 * QREFvis3d(ig,l,iaer) /
402	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
403	& pq(ig,l,igcm_dust_submicron) *
404	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
405	ENDDO
406	ENDIF
407	ENDDO
408
409	c==================================================================
410	CASE("h2o_ice") aerkind ! Water ice crystals
411	c==================================================================
412
413	c 1. Initialization
414	aerosol(1:ngrid,1:nlayer,iaer) = 0.
415	taucloudvis(1:ngrid) = 0.
416	taucloudtes(1:ngrid) = 0.
417	c 2. Opacity calculation
418	DO ig=1, ngrid
419	DO l=1,nlayer
420	aerosol(ig,l,iaer) = max(1E-20,
421	& ( 0.75 * QREFvis3d(ig,l,iaer) /
422	& ( rho_ice * reffrad(ig,l,iaer) ) ) *
423	& pq(ig,l,i_ice) *
424	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
425	& )
426	taucloudvis(ig) = taucloudvis(ig) + aerosol(ig,l,iaer)
427	taucloudtes(ig) = taucloudtes(ig) + aerosol(ig,l,iaer)*
428	& QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer) *
429	& ( 1.E0 - omegaREFir3d(ig,l,iaer) )
430	ENDDO
431	ENDDO
432	c 3. Outputs -- Now done in physiq.F
433	! IF (ngrid.NE.1) THEN
434	! CALL WRITEDIAGFI(ngrid,'tauVIS','tauext VIS refwvl',
435	! & ' ',2,taucloudvis)
436	! CALL WRITEDIAGFI(ngrid,'tauTES','tauabs IR refwvl',
437	! & ' ',2,taucloudtes)
438	! IF (callstats) THEN
439	! CALL wstats(ngrid,'tauVIS','tauext VIS refwvl',
440	! & ' ',2,taucloudvis)
441	! CALL wstats(ngrid,'tauTES','tauabs IR refwvl',
442	! & ' ',2,taucloudtes)
443	! ENDIF
444	! ELSE
445	! CALL writeg1d(ngrid,1,taucloudtes,'tautes','NU')
446	! ENDIF
447	c==================================================================
448	END SELECT aerkind
449	c -----------------------------------
450	ENDDO ! iaer (loop on aerosol kind)
451
452	c -----------------------------------------------------------------
453	c Rescaling each layer to reproduce the choosen (or assimilated)
454	c dust extinction opacity at visible reference wavelength, which
455	c is originally scaled to an equivalent odpref Pa pressure surface.
456	c -----------------------------------------------------------------
457
458	#ifdef DUSTSTORM
459	c -----------------------------------------------------------------
460	! Calculate reference opacity without perturbation
461	c -----------------------------------------------------------------
462	IF (firstcall) THEN
463	DO iaer=1,naerdust
464	DO l=1,nlayer
465	DO ig=1,ngrid
466	tauref(ig) = tauref(ig) +
467	& aerosol(ig,l,iaerdust(iaer))
468	ENDDO
469	ENDDO
470	ENDDO
471	tauref(:) = tauref(:) * odpref / pplev(:,1)
472
473	c--------------------------------------------------
474	c Get parameters of the opacity perturbation
475	c--------------------------------------------------
476	iaer=1 ! just change dust
477
478	write(,) "Add a local storm ?"
479	localstorm=.true. ! default value
480	call getin("localstorm",localstorm)
481	write(,) " localstorm = ",localstorm
482
483	IF (localstorm) THEN
484	WRITE(,) "********************"
485	WRITE(,) "ADDING A LOCAL STORM"
486	WRITE(,) "********************"
487
488	write(,) "ref opacity of local dust storm"
489	taulocref = 4.25 ! default value
490	call getin("taulocref",taulocref)
491	write(,) " taulocref = ",taulocref
492
493	write(,) "target altitude of local storm (km)"
494	ztoploc = 10.0 ! default value
495	call getin("ztoploc",ztoploc)
496	write(,) " ztoploc = ",ztoploc
497
498	write(,) "radius of dust storm (degree)"
499	radloc = 0.5 ! default value
500	call getin("radloc",radloc)
501	write(,) " radloc = ",radloc
502
503	write(,) "center longitude of storm (deg)"
504	lonloc = 25.0 ! default value
505	call getin("lonloc",lonloc)
506	write(,) " lonloc = ",lonloc
507
508	write(,) "center latitude of storm (deg)"
509	latloc = -2.5 ! default value
510	call getin("latloc",latloc)
511	write(,) " latloc = ",latloc
512
513	write(,) "reff storm (mic) 0. for background"
514	reffstorm = 0.0 ! default value
515	call getin("reffstorm",reffstorm)
516	write(,) " reffstorm = ",reffstorm
517
518	!! LOOP: modify opacity
519	DO ig=1,ngrid
520
521	!! distance to the center:
522	ray(ig)=SQRT((lati(ig)180./pi-latloc)*2 +
523	& (long(ig)180./pi -lonloc)*2)
524
525	!! transition factor for storm
526	!! factor is hardcoded -- increase it to steepen
527	yeah = (TANH(2.+(radloc-ray(ig))*10.)+1.)/2.
528
529	!! new opacity field
530	!! -- add an opacity set to taulocref
531	!! -- the additional reference opacity will
532	!! thus be taulocref*odpref/pplev
533	tauuser(ig)=max( tauref(ig) * pplev(ig,1) /odpref ,
534	& taulocref * yeah )
535
536	!! compute l_top
537	DO l=1,nlayer
538	zalt(ig,l) = LOG( pplev(ig,1)/pplev(ig,l) )
539	& / g / 44.01
540	& * 8.31 * 210.
541	IF ( (ztoploc .lt. zalt(ig,l) )
542	& .and. (ztoploc .gt. zalt(ig,l-1)) ) l_top=l-1
543	ENDDO
544
545	!! change reffrad if ever needed
546	IF (reffstorm .gt. 0.) THEN
547	DO l=1,nlayer
548	IF (l .lt. l_top+1) THEN
549	reffrad(ig,l,iaer) = max( reffrad(ig,l,iaer), reffstorm
550	& * 1.e-6 * yeah )
551	ENDIF
552	ENDDO
553	ENDIF
554
555	ENDDO
556	!! END LOOP
557
558	!! compute perturbation in each layer (equation 8 in Spiga et al. JGR 2013)
559	DO ig=1,ngrid
560	int_factor(ig)=0.
561	DO l=1,nlayer
562	IF (l .lt. l_top+1) THEN
563	int_factor(ig) =
564	& int_factor(ig) +
565	& ( 0.75 * QREFvis3d(ig,l,iaer) /
566	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
567	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
568	ENDIF
569	ENDDO
570	DO l=1, nlayer
571	!! Mass mixing ratio perturbation due to local dust storm in each layer
572	more_dust(ig,l,1)=
573	& (tauuser(ig)-(tauref(ig)
574	& * pplev(ig,1) /odpref)) /
575	& int_factor(ig)
576	more_dust(ig,l,2)=
577	& (tauuser(ig)-(tauref(ig) *
578	& pplev(ig,1) /odpref))
579	& / int_factor(ig) *
580	& ((ref_r0/reffrad(ig,l,iaer))**3)
581	& * r3n_q
582	ENDDO
583	ENDDO
584
585	!! quantity of dust for each layer with the addition of the perturbation
586	DO l=1, l_top
587	pq(:,l,igcm_dust_mass)= pq(:,l,igcm_dust_mass)
588	. + more_dust(:,l,1)
589	pq(:,l,igcm_dust_number)= pq(:,l,igcm_dust_number)
590	. + more_dust(:,l,2)
591	ENDDO
592	ENDIF !! IF (localstorm)
593	tauref(:)=0.
594	ENDIF !! IF (firstcall)
595	#endif
596
597	IF (freedust) THEN
598	tauscaling(:) = 1.
599
600	ELSE
601	c Temporary scaling factor
602	taudusttmp(1:ngrid)=0.
603	DO iaer=1,naerdust
604	DO l=1,nlayer
605	DO ig=1,ngrid
606	c Scaling factor
607	taudusttmp(ig) = taudusttmp(ig) +
608	& aerosol(ig,l,iaerdust(iaer))
609	ENDDO
610	ENDDO
611	ENDDO
612
613	c Saved scaling factor
614	DO ig=1,ngrid
615	tauscaling(ig) = tauref(ig) *
616	& pplev(ig,1) / odpref / taudusttmp(ig)
617	ENDDO
618
619	ENDIF ! IF (freedust)
620
621	c Opacity computation
622	DO iaer=1,naerdust
623	DO l=1,nlayer
624	DO ig=1,ngrid
625	aerosol(ig,l,iaerdust(iaer)) = max(1E-20,
626	& aerosol(ig,l,iaerdust(iaer))* tauscaling(ig))
627	ENDDO
628	ENDDO
629	ENDDO
630
631	DO l=1,nlayer
632	DO ig=1,ngrid
633	dsodust(ig,l) = max(1E-20,dsodust(ig,l)* tauscaling(ig))
634	ENDDO
635	ENDDO
636
637	IF (freedust) THEN
638	! tauref has been initialized to 0 before.
639	DO iaer=1,naerdust
640	DO l=1,nlayer
641	DO ig=1,ngrid
642	#ifdef DUSTSTORM
643	!! recalculate opacity because storm perturbation has been added
644	IF (firstcall) THEN
645	aerosol(ig,l,iaer) =
646	& ( 0.75 * QREFvis3d(ig,l,iaer) /
647	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
648	& pq(ig,l,igcm_dust_mass) *
649	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
650	ENDIF
651	#endif
652	tauref(ig) = tauref(ig) +
653	& aerosol(ig,l,iaerdust(iaer))
654	ENDDO
655	ENDDO
656	ENDDO
657	tauref(:) = tauref(:) * odpref / pplev(:,1)
658	ENDIF
659
660	c output for debug
661	c IF (ngrid.NE.1) THEN
662	c CALL WRITEDIAGFI(ngrid,'taudusttmp','virtual tau dust',
663	c & '#',2,taudusttmp)
664	c CALL WRITEDIAGFI(ngrid,'tausca','tauscaling',
665	c & '#',2,tauscaling)
666	c ELSE
667	c CALL WRITEDIAGFI(ngrid,'taudusttmp','virtual tau dust',
668	c & '#',0,taudusttmp)
669	c CALL WRITEDIAGFI(ngrid,'tausca','tauscaling',
670	c & '#',0,tauscaling)
671	c ENDIF
672	c -----------------------------------------------------------------
673	c Column integrated visible optical depth in each point
674	c -----------------------------------------------------------------
675	DO iaer=1,naerkind
676	do l=1,nlayer
677	do ig=1,ngrid
678	tau(ig,iaer) = tau(ig,iaer) + aerosol(ig,l,iaer)
679	end do
680	end do
681	ENDDO
682
683	#ifdef DUSTSTORM
684	IF (firstcall) THEN
685	firstcall=.false.
686	ENDIF
687	#endif
688
689	c -----------------------------------------------------------------
690	c Density scaled opacity and column opacity output
691	c -----------------------------------------------------------------
692	c dsodust(1:ngrid,1:nlayer) = 0.
693	c DO iaer=1,naerdust
694	c DO l=1,nlayer
695	c DO ig=1,ngrid
696	c dsodust(ig,l) = dsodust(ig,l) +
697	c & aerosol(ig,l,iaerdust(iaer)) * g /
698	c & (pplev(ig,l) - pplev(ig,l+1))
699	c ENDDO
700	c ENDDO
701	c IF (ngrid.NE.1) THEN
702	c write(txt2,'(i1.1)') iaer
703	c call WRITEDIAGFI(ngrid,'taudust'//txt2,
704	c & 'Dust col opacity',
705	c & ' ',2,tau(1,iaerdust(iaer)))
706	c IF (callstats) THEN
707	c CALL wstats(ngrid,'taudust'//txt2,
708	c & 'Dust col opacity',
709	c & ' ',2,tau(1,iaerdust(iaer)))
710	c ENDIF
711	c ENDIF
712	c ENDDO
713
714	c IF (ngrid.NE.1) THEN
715	c CALL WRITEDIAGFI(ngrid,'dsodust','tau*g/dp',
716	c & 'm2.kg-1',3,dsodust)
717	c IF (callstats) THEN
718	c CALL wstats(ngrid,'dsodust',
719	c & 'tau*g/dp',
720	c & 'm2.kg-1',3,dsodust)
721	c ENDIF
722	c ELSE
723	c CALL WRITEDIAGFI(ngrid,"dsodust","dsodust","m2.kg-1",1,
724	c & dsodust)
725	c ENDIF ! of IF (ngrid.NE.1)
726	c -----------------------------------------------------------------
727	return
728	end

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

Download in other formats:

Original Format