Context Navigation

aeropacity.F @ 1376

Last change on this file since 1376 was 1376, checked in by aslmd, 10 years ago
amended previous commit: line too long
File size: 26.4 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	REAL zdqnorm(ngrid, nlayer,2)
126	#endif
127
128	c local saved variables
129	c ---------------------
130
131	c Level under which the dust mixing ratio is held constant
132	c when computing the dust opacity in each layer
133	c (this applies when doubleq and active are true)
134	INTEGER, PARAMETER :: cstdustlevel0 = 7
135	INTEGER, SAVE :: cstdustlevel
136
137	LOGICAL,SAVE :: firstcall=.true.
138
139	! indexes of water ice and dust tracers:
140	INTEGER,SAVE :: i_ice=0 ! water ice
141	real,parameter :: odpref=610. ! DOD reference pressure (Pa)
142	CHARACTER(LEN=20) :: txt ! to temporarly store text
143	CHARACTER(LEN=1) :: txt2 ! to temporarly store text
144	! indexes of dust scatterers:
145	INTEGER,SAVE :: naerdust ! number of dust scatterers
146
147	tau(1:ngrid,1:naerkind)=0
148
149	! identify tracers
150
151	IF (firstcall) THEN
152	! identify scatterers that are dust
153	naerdust=0
154	DO iaer=1,naerkind
155	txt=name_iaer(iaer)
156	IF (txt(1:4).eq."dust") THEN
157	naerdust=naerdust+1
158	iaerdust(naerdust)=iaer
159	ENDIF
160	ENDDO
161	! identify tracers which are dust
162	i=0
163	DO iq=1,nq
164	txt=noms(iq)
165	IF (txt(1:4).eq."dust") THEN
166	i=i+1
167	nqdust(i)=iq
168	ENDIF
169	ENDDO
170
171	IF (water.AND.activice) THEN
172	i_ice=igcm_h2o_ice
173	write(,) "aeropacity: i_ice=",i_ice
174	ENDIF
175
176	c typical profile of solsir and (1-w)^(-1):
177	msolsir(1:nlayer,1:naerkind)=0
178	mqextsqabs(1:nlayer,1:naerkind)=0
179	WRITE(,) "Typical profiles of Qext(vis)/Qext(IR)"
180	WRITE(,) " and Qext(IR)/Qabs(IR):"
181	DO iaer = 1, naerkind ! Loop on aerosol kind
182	WRITE(,) "Aerosol # ",iaer
183	DO l=1,nlayer
184	DO ig=1,ngrid
185	msolsir(l,iaer)=msolsir(l,iaer)+
186	& QREFvis3d(ig,l,iaer)/
187	& QREFir3d(ig,l,iaer)
188	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)+
189	& (1.E0-omegaREFir3d(ig,l,iaer))**(-1)
190	ENDDO
191	msolsir(l,iaer)=msolsir(l,iaer)/REAL(ngrid)
192	mqextsqabs(l,iaer)=mqextsqabs(l,iaer)/REAL(ngrid)
193	ENDDO
194	WRITE(,) "solsir: ",msolsir(:,iaer)
195	WRITE(,) "Qext/Qabs(IR): ",mqextsqabs(:,iaer)
196	ENDDO
197
198	! load value of tauvis from callphys.def (if given there,
199	! otherwise default value read from starfi.nc file will be used)
200	call getin("tauvis",tauvis)
201
202	IF (freedust) THEN
203	cstdustlevel = 1
204	ELSE
205	cstdustlevel = cstdustlevel0
206	ENDIF
207
208
209	#ifndef DUSTSTORM
210	firstcall=.false.
211	#endif
212
213	END IF
214
215	c Vertical column optical depth at "odpref" Pa
216	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
217	IF(freedust) THEN
218	tauref(:) = 0. ! tauref is computed after, instead of being forced
219
220	ELSE IF(iaervar.eq.1) THEN
221	do ig=1, ngrid
222	tauref(ig)=max(tauvis,1.e-9) ! tauvis=cste (set in callphys.def
223	! or read in starfi
224	end do
225	ELSE IF (iaervar.eq.2) THEN ! << "Viking" Scenario>>
226
227	tauref(1) = 0.7+.3cos(ls+80.pi/180.) ! like seen by VL1
228	do ig=2,ngrid
229	tauref(ig) = tauref(1)
230	end do
231
232	ELSE IF (iaervar.eq.3) THEN ! << "MGS" scenario >>
233
234	taueq= 0.2 +(0.5-0.2) (cos(0.5(ls-4.363)))**14
235	tauS= 0.1 +(0.5-0.1) (cos(0.5(ls-4.363)))**14
236	tauN = 0.1
237	c if (peri_day.eq.150) then
238	c tauS=0.1
239	c tauN=0.1 +(0.5-0.1) (cos(0.5(ls+pi-4.363)))**14
240	c taueq= 0.2 +(0.5-0.2) (cos(0.5(ls+pi-4.363)))**14
241	c endif
242	do ig=1,ngrid
243	if (lati(ig).ge.0) then
244	! Northern hemisphere
245	tauref(ig)= tauN +
246	& (taueq-tauN)0.5(1+tanh((45-lati(ig)180./pi)6/60))
247	else
248	! Southern hemisphere
249	tauref(ig)= tauS +
250	& (taueq-tauS)0.5(1+tanh((45+lati(ig)180./pi)6/60))
251	endif
252	enddo ! of do ig=1,ngrid
253	ELSE IF (iaervar.eq.5) THEN ! << Escalier Scenario>>
254	c tauref(1) = 0.2
255	c if ((ls.ge.210.pi/180.).and.(ls.le.330.pi/180.))
256	c & tauref(1) = 2.5
257	tauref(1) = 2.5
258	if ((ls.ge.30.pi/180.).and.(ls.le.150.pi/180.))
259	& tauref(1) = .2
260
261	do ig=2,ngrid
262	tauref(ig) = tauref(1)
263	end do
264	ELSE IF ((iaervar.ge.6).and.(iaervar.le.8)) THEN
265	! clim, cold or warm synthetic scenarios
266	call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
267	ELSE IF ((iaervar.ge.24).and.(iaervar.le.31))
268	& THEN ! << MY... dust scenarios >>
269	call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
270	ELSE IF ((iaervar.eq.4).or.
271	& ((iaervar.ge.124).and.(iaervar.le.126))) THEN
272	! "old" TES assimation dust scenario (values at 700Pa in files!)
273	call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
274	ELSE
275	stop 'problem with iaervar in aeropacity.F'
276	ENDIF
277
278	c -----------------------------------------------------------------
279	c Computing the opacity in each layer
280	c -----------------------------------------------------------------
281
282	DO iaer = 1, naerkind ! Loop on aerosol kind
283	c --------------------------------------------
284	aerkind: SELECT CASE (name_iaer(iaer))
285	c==================================================================
286	CASE("dust_conrath") aerkind ! Typical dust profile
287	c==================================================================
288
289	c Altitude of the top of the dust layer
290	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
291	zlsconst=SIN(ls-2.76)
292	if (iddist.eq.1) then
293	do ig=1,ngrid
294	topdust(ig)=topdustref ! constant dust layer top
295	end do
296
297	else if (iddist.eq.2) then ! "Viking" scenario
298	do ig=1,ngrid
299	! altitude of the top of the aerosol layer (km) at Ls=2.76rad:
300	! in the Viking year scenario
301	topdust0(ig)=60. -22.sinlat(ig)*2
302	topdust(ig)=topdust0(ig)+18.*zlsconst
303	end do
304
305	else if(iddist.eq.3) then !"MGS" scenario
306	do ig=1,ngrid
307	topdust(ig)=60.+18.*zlsconst
308	& -(32+18zlsconst)sin(lati(ig))**4
309	& - 8zlsconst(sin(lati(ig)))**5
310	end do
311	endif
312
313	c Optical depth in each layer :
314	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
315	if(iddist.ge.1) then
316
317	expfactor=0.
318	DO l=1,nlayer
319	DO ig=1,ngrid
320	c Typical mixing ratio profile
321	if(pplay(ig,l).gt.odpref
322	$ /(988.**(topdust(ig)/70.))) then
323	zp=(odpref/pplay(ig,l))**(70./topdust(ig))
324	expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
325	else
326	expfactor=1.e-3
327	endif
328	c Vertical scaling function
329	aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) *
330	& expfactor *
331	& QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer)
332	ENDDO
333	ENDDO
334
335	else if(iddist.eq.0) then
336	c old dust vertical distribution function (pollack90)
337	DO l=1,nlayer
338	DO ig=1,ngrid
339	zp=odpref/pplay(ig,l)
340	aerosol(ig,l,1)= tauref(ig)/odpref *
341	s (pplev(ig,l)-pplev(ig,l+1))
342	s max( exp(.03(1.-max(zp,1.))) , 1.E-3 )
343	ENDDO
344	ENDDO
345	end if
346
347	c==================================================================
348	CASE("dust_doubleq") aerkind! Two-moment scheme for dust
349	c (transport of mass and number mixing ratio)
350	c==================================================================
351
352	DO l=1,nlayer
353	IF (l.LE.cstdustlevel) THEN
354	c Opacity in the first levels is held constant to
355	c avoid unrealistic values due to constant lifting:
356	DO ig=1,ngrid
357	aerosol(ig,l,iaer) =
358	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
359	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
360	& pq(ig,cstdustlevel,igcm_dust_mass) *
361	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
362	! DENSITY SCALED OPACITY IN INFRARED:
363	dsodust(ig,l) =
364	& ( 0.75 * QREFir3d(ig,cstdustlevel,iaer) /
365	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
366	& pq(ig,cstdustlevel,igcm_dust_mass)
367	ENDDO
368	ELSE
369	DO ig=1,ngrid
370	aerosol(ig,l,iaer) =
371	& ( 0.75 * QREFvis3d(ig,l,iaer) /
372	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
373	& pq(ig,l,igcm_dust_mass) *
374	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
375	! DENSITY SCALED OPACITY IN INFRARED:
376	dsodust(ig,l) =
377	& ( 0.75 * QREFir3d(ig,l,iaer) /
378	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
379	& pq(ig,l,igcm_dust_mass)
380	ENDDO
381	ENDIF
382	ENDDO
383
384	c==================================================================
385	CASE("dust_submicron") aerkind ! Small dust population
386	c==================================================================
387
388	DO l=1,nlayer
389	IF (l.LE.cstdustlevel) THEN
390	c Opacity in the first levels is held constant to
391	c avoid unrealistic values due to constant lifting:
392	DO ig=1,ngrid
393	aerosol(ig,l,iaer) =
394	& ( 0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
395	& ( rho_dust * reffrad(ig,cstdustlevel,iaer) ) ) *
396	& pq(ig,cstdustlevel,igcm_dust_submicron) *
397	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
398	ENDDO
399	ELSE
400	DO ig=1,ngrid
401	aerosol(ig,l,iaer) =
402	& ( 0.75 * QREFvis3d(ig,l,iaer) /
403	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
404	& pq(ig,l,igcm_dust_submicron) *
405	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
406	ENDDO
407	ENDIF
408	ENDDO
409
410	c==================================================================
411	CASE("h2o_ice") aerkind ! Water ice crystals
412	c==================================================================
413
414	c 1. Initialization
415	aerosol(1:ngrid,1:nlayer,iaer) = 0.
416	taucloudvis(1:ngrid) = 0.
417	taucloudtes(1:ngrid) = 0.
418	c 2. Opacity calculation
419	DO ig=1, ngrid
420	DO l=1,nlayer
421	aerosol(ig,l,iaer) = max(1E-20,
422	& ( 0.75 * QREFvis3d(ig,l,iaer) /
423	& ( rho_ice * reffrad(ig,l,iaer) ) ) *
424	& pq(ig,l,i_ice) *
425	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
426	& )
427	taucloudvis(ig) = taucloudvis(ig) + aerosol(ig,l,iaer)
428	taucloudtes(ig) = taucloudtes(ig) + aerosol(ig,l,iaer)*
429	& QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer) *
430	& ( 1.E0 - omegaREFir3d(ig,l,iaer) )
431	ENDDO
432	ENDDO
433	c 3. Outputs -- Now done in physiq.F
434	! IF (ngrid.NE.1) THEN
435	! CALL WRITEDIAGFI(ngrid,'tauVIS','tauext VIS refwvl',
436	! & ' ',2,taucloudvis)
437	! CALL WRITEDIAGFI(ngrid,'tauTES','tauabs IR refwvl',
438	! & ' ',2,taucloudtes)
439	! IF (callstats) THEN
440	! CALL wstats(ngrid,'tauVIS','tauext VIS refwvl',
441	! & ' ',2,taucloudvis)
442	! CALL wstats(ngrid,'tauTES','tauabs IR refwvl',
443	! & ' ',2,taucloudtes)
444	! ENDIF
445	! ELSE
446	! CALL writeg1d(ngrid,1,taucloudtes,'tautes','NU')
447	! ENDIF
448	c==================================================================
449	END SELECT aerkind
450	c -----------------------------------
451	ENDDO ! iaer (loop on aerosol kind)
452
453	c -----------------------------------------------------------------
454	c Rescaling each layer to reproduce the choosen (or assimilated)
455	c dust extinction opacity at visible reference wavelength, which
456	c is originally scaled to an equivalent odpref Pa pressure surface.
457	c -----------------------------------------------------------------
458
459	#ifdef DUSTSTORM
460	c -----------------------------------------------------------------
461	c the quantity of dust to add at the first time step is calculated
462	c to match a tunable opacity perturbation.
463	c -----------------------------------------------------------------
464	IF (firstcall) THEN
465	DO iaer=1,naerdust
466	DO l=1,nlayer
467	DO ig=1,ngrid
468	tauref(ig) = tauref(ig) +
469	& aerosol(ig,l,iaerdust(iaer))
470	ENDDO
471	ENDDO
472	ENDDO
473	tauref(:) = tauref(:) * odpref / pplev(:,1)
474	c--------------------------------------------------
475	c Parameters of the opacity perturbation
476	c--------------------------------------------------
477
478	iaer=1 !!!! PROVISOIRE !!!!
479
480	write(,) "Add a local storm ?"
481	localstorm=.true. ! default value
482	call getin("localstorm",localstorm)
483	write(,) " localstorm = ",localstorm
484
485	IF (localstorm) THEN
486	WRITE(,) "********************"
487	WRITE(,) "ADDING A LOCAL STORM"
488	WRITE(,) "********************"
489
490	write(,) "ref opacity of local dust storm"
491	taulocref = 4.25 ! default value
492	call getin("taulocref",taulocref)
493	write(,) " taulocref = ",taulocref
494
495	write(,) "target altitude of local storm (km)"
496	ztoploc = 10.0 ! default value
497	call getin("ztoploc",ztoploc)
498	write(,) " ztoploc = ",ztoploc
499
500	write(,) "radius of dust storm (degree)"
501	radloc = 0.5 ! default value
502	call getin("radloc",radloc)
503	write(,) " radloc = ",radloc
504
505	write(,) "center longitude of storm (deg)"
506	lonloc = 25.0 ! default value
507	call getin("lonloc",lonloc)
508	write(,) " lonloc = ",lonloc
509
510	write(,) "center latitude of storm (deg)"
511	latloc = -2.5 ! default value
512	call getin("latloc",latloc)
513	write(,) " latloc = ",latloc
514
515	write(,) "reff storm (mic) 0. for background"
516	reffstorm = 0.0 ! default value
517	call getin("reffstorm",reffstorm)
518	write(,) " reffstorm = ",reffstorm
519
520	DO ig=1,ngrid
521
522
523	c---------------------------------------
524	c distance to the center:
525	c-----------------------------------------
526
527	ray(ig)=SQRT((lati(ig)180./pi-latloc)*2 +
528	& (long(ig)180./pi -lonloc)*2)
529
530
531	!! transition factor for storm
532	!! increase factor ray diff for steepness
533	yeah = (TANH(2.+(radloc-ray(ig))*10.)+1.)/2.
534
535	c-------------------------------------------------
536	c Tau's new map:
537	c------------------------------------------
538
539	tauuser(ig)=max(tauref(ig) * pplev(ig,1) /odpref ,
540	& taulocref * yeah)
541
542	c---------------------------------------------------------
543	c compute l_top
544	c----------------------------------------------------------
545
546	DO l=1,nlayer
547	zalt(ig,l) = LOG( pplev(ig,1)/pplev(ig,l) )
548	& / g / 44.01
549	& * 8.31 * 210.
550	IF ( (ztoploc .lt. zalt(ig,l) )
551	& .and. (ztoploc .gt. zalt(ig,l-1)) ) l_top=l-1
552	ENDDO
553
554	c---------------------------------------------------------
555	c change reffrad if ever needed
556	c----------------------------------------------------------
557
558	IF (reffstorm .gt. 0.) THEN
559	DO l=1,nlayer
560	IF (l .lt. l_top+1) THEN
561	reffrad(ig,l,iaer) = max( reffrad(ig,l,iaer), reffstorm
562	& * 1.e-6 * yeah )
563	ENDIF
564	ENDDO
565	ENDIF
566
567
568	c---------------------------------------------------------------------------
569	ENDDO !! boucle sur ngridmx
570	c---------------------------------------------------------------------------
571	c----------------------------------------------------------------------------
572	c compute int_factor
573	c---------------------------------------------------------------------------
574
575	DO ig=1,ngrid
576	int_factor(ig)=0.
577	DO l=1,nlayer
578	IF (l .lt. l_top+1) THEN
579	int_factor(ig) =
580	& int_factor(ig) +
581	& ( 0.75 * QREFvis3d(ig,l,iaer) /
582	& ( rho_dust * reffrad(ig,l,iaer) ) ) *
583	& ( pplev(ig,l) - pplev(ig,l+1) ) / g
584	ENDIF
585	ENDDO
586	DO l=1, nlayer
587
588	c-------------------------------------------------------------------------
589	c Mass mixing ratio perturbation due to the local dust storm in each
590	c layer
591	c-------------------------------------------------------------------------
592	more_dust(ig,l,1)=
593	& (tauuser(ig)-(tauref(ig)
594	& * pplev(ig,1) /odpref)) /
595	& int_factor(ig)
596	more_dust(ig,l,2)=
597	& (tauuser(ig)-(tauref(ig) *
598	& pplev(ig,1) /odpref))
599	& / int_factor(ig) *
600	& ((ref_r0/reffrad(ig,l,iaer))**3)
601	& * r3n_q
602	ENDDO
603	ENDDO
604
605	c--------------------------------------------------------------------------------------
606	c quantity of dust which is added at the first time step
607	c in dynamical core.
608	c--------------------------------------------------------------------------------------
609
610	DO l=1, l_top
611	zdqnorm(:,l,1) = more_dust(:,l,1)
612	zdqnorm(:,l,2) = more_dust(:,l,2)
613	pq(:,l,igcm_dust_mass)= pq(:,l,igcm_dust_mass)
614	. + zdqnorm(:,l,1)
615	pq(:,l,igcm_dust_number)= pq(:,l,igcm_dust_number)
616	. + zdqnorm(:,l,2)
617	ENDDO
618	ENDIF
619	tauref(:)=0.
620	ENDIF !firstcall
621	#endif
622
623	IF (freedust) THEN
624	tauscaling(:) = 1.
625
626	ELSE
627	c Temporary scaling factor
628	taudusttmp(1:ngrid)=0.
629	DO iaer=1,naerdust
630	DO l=1,nlayer
631	DO ig=1,ngrid
632	c Scaling factor
633	taudusttmp(ig) = taudusttmp(ig) +
634	& aerosol(ig,l,iaerdust(iaer))
635	ENDDO
636	ENDDO
637	ENDDO
638
639	c Saved scaling factor
640	DO ig=1,ngrid
641	tauscaling(ig) = tauref(ig) *
642	& pplev(ig,1) / odpref / taudusttmp(ig)
643	ENDDO
644
645	ENDIF
646
647	c Opacity computation
648	DO iaer=1,naerdust
649	DO l=1,nlayer
650	DO ig=1,ngrid
651	aerosol(ig,l,iaerdust(iaer)) = max(1E-20,
652	& aerosol(ig,l,iaerdust(iaer))* tauscaling(ig))
653	ENDDO
654	ENDDO
655	ENDDO
656
657	DO l=1,nlayer
658	DO ig=1,ngrid
659	dsodust(ig,l) = max(1E-20,dsodust(ig,l)* tauscaling(ig))
660	ENDDO
661	ENDDO
662
663	IF (freedust) THEN
664	! tauref has been initialized to 0 before.
665	DO iaer=1,naerdust
666	DO l=1,nlayer
667	DO ig=1,ngrid
668	tauref(ig) = tauref(ig) +
669	& aerosol(ig,l,iaerdust(iaer))
670	ENDDO
671	ENDDO
672	ENDDO
673	tauref(:) = tauref(:) * odpref / pplev(:,1)
674	ENDIF
675
676	c output for debug
677	c IF (ngrid.NE.1) THEN
678	c CALL WRITEDIAGFI(ngrid,'taudusttmp','virtual tau dust',
679	c & '#',2,taudusttmp)
680	c CALL WRITEDIAGFI(ngrid,'tausca','tauscaling',
681	c & '#',2,tauscaling)
682	c ELSE
683	c CALL WRITEDIAGFI(ngrid,'taudusttmp','virtual tau dust',
684	c & '#',0,taudusttmp)
685	c CALL WRITEDIAGFI(ngrid,'tausca','tauscaling',
686	c & '#',0,tauscaling)
687	c ENDIF
688	c -----------------------------------------------------------------
689	c Column integrated visible optical depth in each point
690	c -----------------------------------------------------------------
691	DO iaer=1,naerkind
692	do l=1,nlayer
693	do ig=1,ngrid
694	tau(ig,iaer) = tau(ig,iaer) + aerosol(ig,l,iaer)
695	end do
696	end do
697	ENDDO
698
699	#ifdef DUSTSTORM
700	IF (firstcall) THEN
701	firstcall=.false.
702	ENDIF
703	#endif
704
705	c -----------------------------------------------------------------
706	c Density scaled opacity and column opacity output
707	c -----------------------------------------------------------------
708	c dsodust(1:ngrid,1:nlayer) = 0.
709	c DO iaer=1,naerdust
710	c DO l=1,nlayer
711	c DO ig=1,ngrid
712	c dsodust(ig,l) = dsodust(ig,l) +
713	c & aerosol(ig,l,iaerdust(iaer)) * g /
714	c & (pplev(ig,l) - pplev(ig,l+1))
715	c ENDDO
716	c ENDDO
717	c IF (ngrid.NE.1) THEN
718	c write(txt2,'(i1.1)') iaer
719	c call WRITEDIAGFI(ngrid,'taudust'//txt2,
720	c & 'Dust col opacity',
721	c & ' ',2,tau(1,iaerdust(iaer)))
722	c IF (callstats) THEN
723	c CALL wstats(ngrid,'taudust'//txt2,
724	c & 'Dust col opacity',
725	c & ' ',2,tau(1,iaerdust(iaer)))
726	c ENDIF
727	c ENDIF
728	c ENDDO
729
730	c IF (ngrid.NE.1) THEN
731	c CALL WRITEDIAGFI(ngrid,'dsodust','tau*g/dp',
732	c & 'm2.kg-1',3,dsodust)
733	c IF (callstats) THEN
734	c CALL wstats(ngrid,'dsodust',
735	c & 'tau*g/dp',
736	c & 'm2.kg-1',3,dsodust)
737	c ENDIF
738	c ELSE
739	c CALL WRITEDIAGFI(ngrid,"dsodust","dsodust","m2.kg-1",1,
740	c & dsodust)
741	c ENDIF ! of IF (ngrid.NE.1)
742	c -----------------------------------------------------------------
743	return
744	end

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

Download in other formats:

Original Format