Context Navigation

aeropacity_mod.F @ 1944

Last change on this file since 1944 was 1861, checked in by emillour, 7 years ago
Mars GCM: Add possibility to load a dust MY33 scenario. EM
File size: 27.4 KB

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

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

Download in other formats:

Original Format